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THESIS DISCLAIMER 


The reader is cautioned that computer programs developed in this research may 
not have been exercised for all cases of interest. While every effort has been made, 
vvithin the time available, to ensure that the programs are free of computational and logic 
errors, they cannot be considered validated. Any application of these programs without 
additional verification is at the risk of the user. 
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EXECUTIVE SUMMARY 


The passive location of both hostile and friendly electromagnetic emitters has been 
an important capability for the war-fighter, for law enforcement, and in search and rescue 
operations. The characteristics of an emitted signal once received by several sensors in an 
operating environment can be exploited to passively locate an emitter. Two of the 
characteristics that can be measured by separate receivers of such a signal are: (1) the 
differences in Doppler shifts (if any) in the signal’s frequency and (2) the difference in the 
time that the signal arrived at the sensors. These measurements are known as Frequency 
Difference of Arrival (FDOA) and Time Difference of Arrival (TDOA). In the presence of 
moving sensors, the sensors’ receivers measure the frequency of arrival (FOA) of the 
signal. The difference of two FOA’s provides an FDOA which yields a surface called an 
isodop, in this case meaning equally differenced Doppler data, which contains the locus of 
all points on which the emitter could lie given the Doppler information. Similarly, the 
difference of a single receiver’s time of arrival (TO A) with another TO A yields a TDOA 
which describes a surface known as an isochron, in this case meaning equally differenced 
time, which contains the locus of all points on which the emitter could lie given the TDOA 
information. 

In the past, the measurement of such quantities in the tactical environment has 
been difficult due to their sensitivities to timing errors between any two sensors. Further, 
when trying to use these measurements for the geolocation of an emitter, errors in the 
sensor’s own location and velocity measurements compound with the signal processing 
errors to provide unreliable emitter fixes. Historically, the Navy has relied upon long 
range, shore-based High Frequency Direction Finding (HFDF) networks. These networks 
use large antennas to measure the angle of arrival of a signal and cross these bearings over 
long distances to triangularize the emitter’s location. Given the recent losses of Navy 
overseas assets and that more and more emitters in the tactical environment operate in the 
Very High Frequency (VHF) and Ultra High Frequency (UHF) range, the ability to 
perform emitter geolocation is diminishing. As a result, the Naval Security Group Support 
Activity (NSGSA) contracted with the Applied Research Lab at the University of Texas at 
Austin (ARL:UT) to develop an affordable, low-risk TDOA geolocation system using 
commercial off the shelf (COTS) and government off the shelf (GOTS) technologies. 
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ARL;UT responded by the development of the Carry-On Multi-Platform Global 
Positioning System Assisted TDOA System. 

When using TDOA techniques, errors in the emitter geolocation problem is 
essentially related to three components: the location of the observer relative to the emitter 
(the “geometry”), signal timing, and observer position measurement. The emergence of 
Global Positioning System (GPS) technology offers capabilities which can greatly enhance 
TDOA applications by reducing measurement error components. Specifically, GPS allows 
for dramatic reductions in the observer position error and sample timing between 
observers. In addition, the advantages of using such commercially available and 
government "off-the-shelf technology allows the system to be placed on a variety of 
platforms and, thus, could make the GPS-Assisted TDOA Geolocation System a highly 
desirable Multi-Platform tool in the "Joint Warfare" environment. 

The measurement of the FDOA is done simultaneously when measuring the TDOA 
in the formulation used by the ARL:UT Test System. While the initial test and 
development of the ARL:UT system focused on the use of only TDOA measurements for 
geolocation, the combined use of TDOA and FDOA measurements for geolocation has 
become possible due to the GPS technology in the ARL:UT system. By providing more 
information to the geolocation solution, combined TDOA and FDOA measurements can 
expand the size of the state of the observer—providing velocity information as well as the 
emitter’s location. In addition to providing velocity information, it has been postulated 
that when poor geometry exists between observers, FDOA can more tightly resolve the 
solution of the geolocation fix than TDOA alone. Further, when the emitter is believed to 
have zero velocity, FDOA can provide the emitter’s three dimensional location using only 
three observers instead of four by fixing the solution’s velocity state to zero. 

Regardless of which method or algorithm is used, the measurement of the TDOA 
and the FDOA is a statistical process, built upon standard estimation assumptions. As 
such, any measurement of either component can be viewed as a random variable related to 
the signal, the environment, and the accuracy and precision of the whole system that 
performs the measurement. Through the use of a Java-based, Stochastic TDOA and 
FDOA Simulation (JBSTAFSim), this thesis explores the impact on the geolocation of an 
emitter in the presence of such measurement errors. Given an estimate of the distribution 
of these measurement errors, JBSTAFSim can simulate the impact on a TDOA or mixed 
TDOA/FDOA geolocation solution within a given “Joint Warfare” scenario of observers. 
In particular, JBSTAFSim allows the user to: adjust the magnitude and correlation of 
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measurement error variance; adjust the believed measurement error variance used by the 
geolocator; dictate the network geometry and types of observers; and elect to choose a 
batch processed geolocator which incorporates previous geolocation information into the 
current solution to solve for the emitter’s state. 

Through simulation, this thesis shows that the mixed T/FDOA solution 
stochastically dominates the TDOA solution in two-dimensional fix scenarios. While no 
stochastic dominance of either solution is shown in the three-dimensional fix case, the 
author shows that the accuracy of the three-dimensional problem is related to the sensor 
network geometry. Further, if a three-dimensional fix is required, it is shown that robust 
sensors like satellites should be used to improve the sensor network to target geometry. 
Finally, since the location of the target relative to the sensors is generally not known, the 
author demonstrates the need to define a method to rate the sensor network geometry 
given possible target locations. 

These results can be viewed fi'om three perspectives. The first view is from the 
perspective of the designer of the system who wishes to understand the sensitivities of 
measurement errors and wishes to improve upon these measurements by the system or 
take them into account in the geolocation process. The next view is fi'om the warrior, or 
user of the system, who would be employing this system to localize a hostile or friendly 
emitter. The warrior is concerned with the quality of the information provided by the 
system. As such, the warrior needs to know how well the system can perform and what 
he or she can do with respect to the sensor network to improve the fix of the target. 
Finally, the last perspective is fi'om that of the activities which fund the designer to create 
such systems for the warrior. This activity is concerned with system performance for the 
amount invested. Given several functional areas that system designers could pursue to 
develop or improve, the simulation can rate the impact of this additional or improved 
information on the system’s performance and describe how this impacts the warfare 
commander’s utility for such information. 
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I. INTRODUCTION 


A. BACKGROUND 

The passive location of both hostile and friendly electromagnetic emitters has been 
an important capability for the war-fighter, for law enforcement, and in search and rescue 
operations. The characteristics of an emitted signal once received by several sensors in an 
operating environment can be exploited to passively locate an emitter. Two of the 
characteristics that can be measured by separate receivers of such a signal are; (1) the 
differences in Doppler shifts (if any) in the signal’s frequency and (2) the difference in the 
time that the signal arrived at the sensors. These measurements are known as Frequency 
Difference of Arrival (FDOA) and Time Difference of Arrival (TDOA). In the presence of 
moving sensors, the sensors’ receivers measure the frequency of arrival (FOA) of the 
signal. The difference of two FOA’s provides an FDOA which yields a surface called an 
isodop, in this case meaning equally differenced Doppler data, which contains the locus of 
all points on which the emitter could lie given the Doppler information (see Figure 1-1). 
Similarly, the difference of a single receiver’s time of arrival (TOA) with another TOA 
yields a TDOA which describes a surface known as an isochron, in this case meaning 
equally differenced time, which contains the locus of all points on which the emitter could 
lie given the TDOA information (see Figure 1-2). 
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Figure 1-2 TDOA’s of Three Observers (From Ref 3) 


In the past, the measurement of such quantities in the tactical environment has 
been difficult due to their sensitivities to timing errors between any two sensors. Further, 
when trying to use these measurements for the geolocation of an emitter, errors in the 
sensor’s own location and velocity measurements compound with the signal processing 
errors to provide unreliable emitter fixes. Historically, the Navy has relied upon long 
range, shore-based High Frequency Direction Finding (HFDF) networks. These networks 
use large antennas to measure the angle of arrival of a signal and cross these bearings over 
long distances to triangularize the emitter’s location. Given the recent losses of Navy 
overseas assets and that more and more emitters in the tactical environment operate in the 
Very High Frequency (VHF) and Ultra High Frequency (UHF) range, the ability to 
perform emitter geolocation is diminishing. As a result, “the Naval Security Group 
Support Activity (NSGSA) contracted with the Applied Research Lab at the University of 
Texas at Austin (ARL;UT) to develop an affordable, low-risk TDOA geolocation system 
using commercial off the shelf (COTS) and government off the shelf (GOTS) 
technologies”[Ref 1]. ARL;UT responded by the development of the Carry-On Multi- 
Platform Global Positioning System Assisted TDOA System. 

When using TDOA techniques, errors in the emitter geolocation problem is 
essentially related to three components; the location of the observer relative to the emitter 
(the “geometry”), signal timing, and observer position measurement. The emergence of 
Global Positioning System (GPS) technology offers capabilities which can greatly enhance 
TDOA applications by reducing measurement error components. Specifically, GPS allows 
for dramatic reductions in the observer position error and sample timing between 
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observers. In addition, the advantages of using such commercially available and 
government "off-the-shelf technology allows the system to be placed on a variety of 
platforms and, thus, could make the GPS-Assisted TDOA Geolocation System a highly 
desirable Multi-Platform tool in the "Joint Warfare" environment as illustrated in Figure 
1-3. 



Figure 1-3 GPS Assisted Geolocation System in Joint Warfare (From Ref 2) 

Through the use of GPS technology and differential techniques combined with Rubidium 
standard oscillators, the total theoretical signal timing error is approximately 25 nano¬ 
seconds which, when multiplied by the speed of light, converts to 7.5 meters in 
distance. [Ref 1] 

The measurement of the FDOA is done simultaneously when measuring the TDOA 
in the formulation used by the ARLiUT Test System While the initial test and 
development of the ARL:UT system focused on the use of only TDOA measurements for 
geolocation, the combined use of TDOA and FDOA measurements for geolocation has 
become possible due to the GPS technology in the ARL:UT system. By providing more 
information to the geolocation solution, combined TDOA and FDOA measurements can 
expand the size of the state of the observer—providing velocity information as well as the 
emitter’s location. In addition to providing velocity information, it has been postulated 
that when poor geometry exists between observers, FDOA can more tightly resolve the 
solution of the geolocation fix than TDOA alone. Further, when the emitter is believed to 
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have zero velocity, FDOA can provide the emitter’s three dimensional location using only 
three observers instead of four by fixing the solution’s velocity state to zero. 


B. PURPOSE 

Regardless what method or algorithm used, the measurement of the TDOA and the 
FDOA is a statistical process, built upon standard estimation assumptions. As such, any 
measurement of either component can be viewed as a random variable related to the 
signal, the signal’s environment, and the accuracy and precision of the whole system that 
performs the measurement. Through the use of a Java-based, Stochastic TDOA and 
FDOA Simulation (JBSTAFSim), the author explores the impact on the geolocation of an 
emitter in the presence of such measurement noise. Given an estimate of the variance of 
these measurement errors, JBSTAFSim can simulate the impact on a TDOA or mixed 
TDOA/FDOA geolocation solution within a given “Joint Warfare” scenario of observers. 
In particular, JBSTAFSim allows the user to: adjust the magnitude and correlation of 
measurement error variance; adjust the believed measurement error variance used by the 
geolocator; dictate the network geometry and types of observers; and elect to choose a 
batch processed geolocator which incorporates previous geolocation information into the 
current solution to solve for the emitter’s state. 

These results can be viewed from three perspectives. The first view is from the 
perspective of the designer of the system and wishes to understand the sensitivities of 
measurement errors and wishes to improve upon these measurements by the system or 
take them into account in the geolocation process. The next view is from the warrior, or 
user of the system, who would be employing this system to localize a hostile or friendly 
emitter. The warrior is concerned with the quality of the information provided by the 
system. As such, the warrior needs to know how well the system can perform and what 
he or she can do with respect to the sensor network to improve the fix of the target. 
Finally, the last perspective is from that of the activities which fund the designer to create 
such systems for the warrior. This activity is concerned with system performance for the 
amount invested. Given several functional areas that system designers could pursue to 
develop or improve, the simulation can rate the impact of this additional or improved 
information on the system’s performance and describe how this impacts the warfare 
commander’s utility for such information 
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C. ORGANIZATION 


The next chapters discuss briefly: the configuration of the ARL:UT test system, 
how the TDOA and FDOA measurements are made, and how these measurements enter 
into the geolocation process. Some of the detailed theory and methods that JBSTAFSim 
uses can be found in Appendixes A-C. Chapters V and VI provide an overview of 
JBSTAFSim and a description of the simulated environment. Chapter VII presents the 
results of the simulation of TDOA and mixed TDOA/FDOA geolocation methods. 
Chapter VIII offers recommendations based on the results presented in the previous 
chapter. 
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n. ARL:UT SYSTEM OVERVIEW 


The ARL:UT system is capable of receiving a frequency in the HF, VHF, or UHF 
frequency bands. Only one voice grade channel compatible with existing Navy hardware 
for communication between sensors is required. Further, the system is designed to be 
compatible with the Unified Build (UB) environment and interface with and display its 
results graphically on the Joint Maritime Command and Information System (JMCIS). 
Each observer updates its location with a GPS receiver. The receivers have a one pulse- 
per-second (pps) output which is used to trigger sampling of the incoming signal of 
interest. The accuracy of the GPS clocks drives each receivers’ output of the one pps 
pulse to be within 20 nano seconds of each other. The sampler at each observer, then, is 
phase-locked to the one pps pulse to within five nano seconds making the samples 
between observers to be within 25 nano seconds [Ref 20]. As depicted in Figure 2-1, the 
system network consists of one master sensor and a few slave sensors. In the evenly 



Figure 2-1 System Overview (From Ref 20) 


determined case, three sensors are required for a two-dimensional geolocation and four for 
a three-dimensional fix. Once given the frequency of interest from traditional detection 
means, the master tunes its receiver to that frequency and orders the slaves to tune to the 
same frequency. Each observer stores the incoming signal as buffered data. Since only 
one communications channel is required to be available, the ARL:UT system uses a 
distributed processing system to make the most efficient use of the single communications 
link. The master determines a 400 milli second portion that contains the most energy and 
computes the Fast Fourier Transform (FFT) of this sample. This FFT is then sent to the 
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slaves with a GPS time stamp. The slaves search their buffers for the 400 milli second 
sample that corresponds to the master’s time stamped FFT. The slaves then FFT their 
own data and correlate these samples using a Complex Ambiguity Function (CAF) to 
determine the TDOA. The TDOA, FDOA, GPS position of the observer, velocity of the 
observer, and timestamp are returned to the master. The master collects the information 
and sends the package to a standard Navy TAC3/4 workstation which calculates the 
geolocation. A detailed reference of the preliminary hardware configuration and 
specifications of the ARL:UT system can be found in Reference 20. 

Like any system that depends on the Global Positioning System, this system is very 
vulnerable if GPS were lost or significantly degraded. While relatively reliable sensor 
information such as position and velocity might be available fi'om inertial systems on the 
sensor’s platform, the need for a reliable, common time standard between sensors so that 
data can be buffered and retrieved at common times would still be required. While the 
vulnerabilities of GPS is an issue well beyond the scope of this discussion, it must be 
considered when weighing the value, limitations, and capabilities of such a system. 
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m. TDOA AND FDOA MEASUREMENTS 


A. BACKGROUND AND PURPOSE 

The main purpose of this chapter is to describe how the TDOA and FDOA 
measurements are made and the assumptions behind those measurements. Much of the 
next section is a basic summary of several signal processing theories and methods. If the 
reader is not interested in these details, the first three sections of this chapter may be 
omitted with the following in mind. The main importance of this chapter is to present the 
measurement of TDOA 's and FDOA’s from a stochastic RF reception as a statistical 
process and to get an understanding of the assumptions required to formulate these 
processes. The measurement errors in any TDOA or FDOA are multivariate functions of 
many parameters—^many of which are random variables themselves. Measurements are 
further confounded by the sensitivity and accuracy of the method employed to determine 
the TDOA or FDOA measurement. 

JBSTAFSim does not assume any one particular method to measure the TDOA or 
FDOA; nor does JBSTAFSim assume any properties with respect to the signal, the noise 
in the signal, or the noise in the environment—the random variables which have been 
shown in the works hsted below to function as parameters which describe the 
measurement error estimate. Rather, JBSTAFSim allows the user to set the measurement 
error variance estimate and correlation of the measurements. Since these measurements 
are modeled statistically, it is reasonable to assume that regardless of the method used to 
obtain them or the many parameters which might describe them, the error in their 
measurement can also be described statistically. By doing this, JBSTAFSim can simulate 
these measurements and indicate their effect on the geolocation process with minimal 
information provided with respect to; the system performing the measurements, the signal 
of interest, and the many random variables describing the received signal and its 
environment. Viewed this way, the model maintains integrity while also allowing it to 
serve in the most general sense. 

The following discussion is based on Reference 1, an NPS thesis by LT David A. 
Streight who applies Cyclostationary techniques to measure the TDOA of an emitter. In 
describing the advantages of this particular method LT Streight provides an excellent 
overview of the basic statistical properties and assumptions of the signal involved in 
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determining the TDOA and FDOA measurement using traditional methods. More 
information about these process can be found in detail in References 4, 5, and 6. 


B. BASIC MODEL AND THE GENERALIZED CROSS CORRELATION 

Most TDOA models assume a stationary signal which is received by two sensors in 
the presence of white Gaussian noise. These signals have been mathematically modeled as 
a function in the time domain: 


x (/) = 5 (/)+ /?,(/) 

2-1 

y{l)= A-s{l-D) + n^{t) 

2-2 


where 5(/) is assumed to be statistically independent and uncorrelated with independent 
noises «,(/) and n^{t). A is the complex valued relative magnitude and phase mismatch 

between the receivers; D represents the TDOA between the signals. [Ref 1] The 
autocorrelation and cross-correlation functions are given by: 

Ki^) = Ki^) + KM) 2-3 

i?^.(r) = Up/?,(r) + 7^_(r) 2-4 


^>.(0 = A-RXt-D) + (r) 


2-5 


and the spectral density functions are: 


5.(/)=5,(/) + S,(/) 

2-6 

5,(/) = Mr5,(/) + .J,(/) 

2-7 


2-8 
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Using these relationships and assumptions, the parameter D, the TDOA, can be 

obtained by using correlation techniques. The most basic technique is the Generalized 
Cross Correlation method. First, note that (by 2-5) peaks when t = D. Next 

assuming that / of the received signal is Avithin a finite Bandwidth B around the carrier 
frequency the ratio of *S'^(/) to 


0 


B B 

/o'y-Z-Zo+Y 
otherwise 


2-9 


The inverse Fourier transform of this ratio is equal to; 


^o(^) 


/o+5/2 ^ (A 

/o-B/2 ) 


which peaks at r = Z) and provides the TDOA measurement. [Ref 1] 
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C. COMPLEX AMBIGUITY FUNCTION 

The Generalized Cross Correlation method as described above works well in the 
case of static transmitters and receivers. However, as shown in Reference 5, in order to 
determine the TDOA when there exists a Doppler shift, the spectrum of one of the signals 
must be translated in fi'equency equal to the FDOA measured between the observers. This 
can be shown by rewriting the generalized model above so that: 


x{t) = s{t) + n,{t) 

2-11 

>^(0 = A-s(t- (0 

2-12 


where is the Doppler difference measured between the two observers. This requires 
the maximization over the two parameters D and which leads to: 
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2-13 


= J 4.‘)y'(.‘ - D)e-'‘^‘'cll 
0 

Notice that at /^ = 0, the CAF reduces to the Generalized Cross Correlation method. In 
Reference 5, Stein shows, given these assumptions, the variance of the estimates for the 
TDOA and FDOA can be related to the signal integration time, the signal bandwidth, and 
the noise bandwidth: 


1 1 

~ p .Jbt^ 


2-14 


1 1 
~ T, .jmy 


2-15 


where: 


B= the noise bandwidth at the input of both receivers. 


A 

P=2k 



\rw,(s)df 


!>K(/)d/ 

_ 


2-16 


with IV ^(/) = the spectral density of the signal as shaped by the receiver. 


=rms. integration time and. 


2 

7 


2 

2 


1 1 

— + — + 

7 . 72 


1 

7x72_ 


2-17 


with — = the signal to noise ratio (SNR) for each receiver [Ref 1] 

X / 
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With a little sensitivity analysis of the above equations, it can be seen that the driving 
component in the standard deviation of the TDOA measurement is SNR. In an idealized 
case, below lOdB SNR and with a 400 ms integration time of the ARL:UT system, the 
standard deviation of a single TDOA measurement exceeds 100 meters. “With three 
measurements to provide a fix, all within this idealized case, the geolocation would have a 
standard deviation of 173 meters” [Ref 1], a large best-case error. 


D, CYCLOSTATIONARY DETERMINATION 

LT Streight’s thesis applies Cyclostationary techniques to measure the TDOA. 
This is a robust approach which is able to perform well despite both low SNR and the 
presence signals which may or may not intentionally seek to jam the observers’ receivers. 
By taking advantage of the cyclic—time periodic—nature of most signals, a three 
dimensional correlation in the fi'equency and cyclic fi'equency domains can be performed. 
Thus, where the signal of interest and noise may have appeared overlapping in the 
traditional models, they can still be separated by a cyclic autocorrelation function which 
can be thought of as a special autocorrelation function which produces spectral lines at 
frequencies indicative of the signal’s periodic nature. In the frequency domain, the cyclic 
autocorrelation function is replaced by the spectral correlation function. The Spectral 
Correlation Function (SCF) is a correlation in the frequency domain which peaks at the 
nominal spectral fi'equency and cyclic frequency in what is called a bi-ffequency plane. 
See Figure 2-1. 
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Figure 3-1 Signals plotted in the Bi-frequency Plane (From Ref. 1) 


This method allows for dramatically better results in the presence of greater 
interference and, as seen above, can provide a classification of the modulation type. The 
model of the signal at any pair of receivers remains as in 2-1 and 2-2, 

x{i) = s(i) + n,{t) 
y{i)= As{t-D) + n^{t), 

and assumes the same statistical assumptions with respect to the properties of the signal 
and noise except that now the noise may now contain co-channel interference which could 
spectrally mask the signal. As before, auto and cross-correlation functions are used in 
conjunction with the spectral density functions to measure the TDOA; however, now 
these correlated values are aligned by the cyclic frequency characteristic of the desired 
signal. While highly accurate in the presence of interference, this method is 
computationally expensive. LT Streight employs a computationally efficient 
cyclostaionary TDOA algorithm, SPECCOA, “which may not be as be as accurate as 
others but is the best choice for the efficiency required in a tactical system.”[Ref 1] While 
not implemented by the ARL:UT test system, this method could be implemented as the 
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algorithm of choice in the next system or in future systems. It is presented here in detail, 
along with the GCC and CAF, to point out the multitude of processes, parameters, and 
assumptions required to estimate the TDOA or FDOA measurement error variance. 
Regardless of the method used to obtain them, the TDOA and FDOA measurement errors 
can be represented as a stochastic process. 


E. IMPLICATIONS OF SIGNAL PROCESSING TO GEOLOCATION 

Throughout this chapter, discussion has focused on the signal models and their 
assumptions. As noted previously Stein, in Reference 5, gives the TDOA and FDOA 
estimate measurement errors based on a few parameters and several assumed conditions. 
Dr. Petre Rusu and Dr. Lisa Giulianelli of ARL:UT, in References 7 and 8, follow the 
work of Azaria and Hertz, in Reference 4, to find the solution for the correlation between 
TDOA measurements as a function of SNR, observation interval length, and noise and 
signal bandwidths. Unfortunately, many of these parameters that go into these 
formulations are in fact stochastic themselves. Further, in Reference 7, the enabling 
assumptions with regard to the estimation of the measurement error correlation begin with 
those of equations 2-1 and 2-2, namely, that the signal received at any observer can be 
described as the summation of two parts—a deterministic signal plus Gaussian noise. The 
noises are considered “uncorrelated with the signal and with each other across 
observers”[Ref 7]. While these assumptions like these are critical for the necessary 
mathematical formulations like those in sections B through D above and in works such as 
Reference 7, in an operating environment, these models will fail to represent the one 
stochastic, and most important, variable present at the observers’ receivers—^the actual, 
“real world” signal itself Viewing each signal received at each of the observers as a 
stochastic process, the TDOA and FDOA measurement process should not be thought of 
as that which can parameterize the measurement errors; rather, it is the randomness 
nascent in each observer’s received signal that introduces the measurement errors. 

Since quantifying the impacts of the estimate covariance and correlation of TDOA 
and FDOA measurement errors on the geolocation process is the goal of these analyses, 
this discussion might seem to be disputing a fine point. However, it demonstrates the 
evaluative value of JBSTAFSim as a tool for assessing the impact of any assumption or 
model of the TDOA and FDOA measurement errors on the geolocation process. 
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JBSTAFSim, then, functions as an “If—Then” tool. The “If—Then” portion of the 
simulation is remarkably simple: given a network and allocation of observers, if 
measurement errors are introduced with a given probabilistic description; then, the 
expected Geolocation results are as follows. JBSTAFSim, then, exploits one of the 
greatest advantages of simulation; rather than redefine our mathematical models or 
simulate only the equations of our models, let us instead violate our own assumptions, 
stochastically, and examine the results within an empirical simulation space. The results of 
such an investigation are useful for anyone who: is designing a system such as this and 
wishes to evaluate or improve its effectiveness; using a system such as this for tactical 
decisions and wishes to know its capabilities and limitations; and for anyone sponsoring 
the development of such a system and wishes to know how added capabilities or 
improvements might improve the warrior’s use of such a system. 
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rv. EMITTER GEOLOCATION 


A. BACKGROUND 

The TDOA and FDOA measurements computed by the signal processing methods 
described in the previous chapter form the basis for the measurement data required to 
solve for the emitter’s location. These quantities measured in seconds and in Hertz can be 
represented formally in terms of length (meters) and velocity (meters/second) by the 
following algebraic equations: 


rD(94 =c(/, -/J 

+[y.-ysy O'. 


4-1 


FDOA, = -fs 


1 

1 

1 

1 

1 

1 

j 



y 


4-2 


Where the subscripts / and j represent the observer pair forming the differences and s 
denotes the signal’s source. An observer’s three dimensional location in space is defined 
by (x,y,z), and r and v represent the radial and velocity vectors of the target or 
observers. Therefore, the TDOA and FDOA can be thought of as spatial differences and 
relative velocity differences, rather than in time and frequency units 

Much has been written and theorized about the best way to solve for the state of 
an emitter given TDOA information. While the TDOA measurement is non-linear in the 
emitter’s state, solutions can be formulated to directly solve for the location of the target 
as linear models which satisfy non-linear constraints. These solutions usually involve 
squaring the TO A equation which results in two locations which satisfy the TDOA data. 
References 9, 10, and 12 all show different, but mathematically equivalent algebraic 
methods to directly solve for the state of an emitter. In particular, in Reference 10, Dr. 
Petre Rusu, of ARL:UT, goes through great work to linearly propagate the estimate 
errors through his solution to provide an uncertainty estimate for his solution to the 
emitter’s state. In Reference 11, Dr. Rusu further proves, given consistent weighting 
matrices, the theoretical equivalence of the TOA direct solution and least squares 
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formulation. References 9 and 13 are specifically geared to the algebraic solution of GPS 
equations; however, their techniques can be applied to any TDOA geolocation problem. 

In JBSTAFSim, the method of Reference 13, by Chaffee and Abel, is used to 
provide an initial starting position for the linearized least squares approach. This method 
is very elegant and numerically robust. Further, it provides a useful understanding of the 
solution ambiguity by not immediately squaring the TO A equation. Thus, this method was 
chosen for its numerical superiority and for its elegant resolution of the ambiguity of 
emitter location. The linearized least squares solution uses this direct solution as its 
starting point from which it seeks a solution to minimize the sum of squared errors. After 
a fix is computed, the solver then uses the previous solution as the starting position for the 
subsequent solutions until a request is made for another direct solution to update the 
current solution. The periodicity of these updates may be set by the user. In an actual 
geo-location system, direct solution updates might be requested as a function of many 
complex variables and parameters, for instance: quality of the previous fix, the 
dimensionless shock, time since the last fix, by request of the system user, or quality of the 
received signal. JBSTAFSim seeks to minimize the sum of squared errors by 
implementing a Square Root Information Filter (SRIF). The SRIF was chosen for its 
numerical superiority and its ability easily include a priori statistics and function as a 
numerically stable extended Kalman filter. The SRJF in JBSTAFSim can perform either as 
a single epoch state solver or as an extended Kalman Filter, “batch”, processor. 
JBSTAFSim’s SRIF does not, however, incorporate a movement model that considers a 
target which moves with time. While a movement model would no doubt provide the 
most robust approach, given limited communications between observers, the complexity 
of such a model, the non-linearity of the measurements, and the purpose of the system to 
act as a locator and not as a fire control system, the batch process mode without a 
movement model should be considered adequate for those targets with zero, or near zero, 
velocities. The development of movement models for particular classes of targets would 
be a large undertaking, however, their design could be facilitated through empirical 
simulation by extending the JBSTAFSim solver. 

Finally, and as will be discussed in the next chapter, JBSTAFSim’s simulation 
results should not be viewed as a measure of the ARL;UT system. In particular, 
ARL:UT’s system is currently in the test and development phase in which its Engineers 
and Scientists are considering a multitude of design issues like the ones presented in this 
thesis. Therefore, JBSTAFSim’s results are not necessarily identical with that of the 
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B. DlKECl GEOLOCATlON SOLUTION 


This niethod was introduced in Reference 13 in order to provide a description of 
the exact solution for the GPS pseudorange equations. Tliis niethod provides a greater 
understanding of the geometiy involved in the TDOA equations and can be applied in a 
nunierically stable fasliion. Dr. Toirnan of ARL:UT provides a simple explaiiation of 

^^iiiuicc ailu /^Dci njiiiiuiauuii ui jvciciciicc ih. ims scCuuii piuviucS a quiCR 5>uiiuiicuy 
of this direct solution method based on both of these references. 


Following the notation and example of Dr. Tolman’s notes in Reference 14, let the 
emitter be located at position F and define an observer to be located at position r.. The 

range between one observer and the emitter is defined by: 


— Iir ^ f II 


4-3 


The TOA equation is defined as: 

TOA, = ct. = - F jj + ct^ =R,+ c/o 4-4 

where 4 is the unknown time of transmission of the signal. Tliis equation holds for any 
change in the origin of time. The TDOA is simply the difference of two TOA’s; 

-«/) 


Algoritlirns in References 9 and 10 begin by squarhig the TOA equation above. Chaffee 
and Abel do not. Instead, label one receiver I, the master. Let ail other receivers be 

_Tx-L : _i ay 1 

UCllOieU VVlLll / 1,... iV'i — 1 . 


c/, - ctj - [ct, - aj) = f, -rj-{r-Fj )j} 


4-6 
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Further, define; 


A, =cr, -ctj 


-n 


r'=r 



This makes the TO A equation; 


A. + r =r’-rii 


Square this equation and note llr iP cancels; 


Irf +2r'*F'= a! +1W'\ 


Stack this equation for M-1 differences yields; 
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A,' 
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4-7 


4-8 


4-9 


For continuity, Chaflfee and Abel describe this set of equations as; 


/r^.i 


U 

\ 1/ 

T5;l 

r5,i\ 

li 

,Z 

JML, 

r^\ 

,k] 

IaJ) 


Where (v,m') represents the Lorentz inner product on 9?“* defined as; if v = (x',a) and 

w = (j',c) are vectors, then {v,w) = x‘y-ac. In addition, also note Z 

Sj =r^', and n = r' [Ref 13] 

Note that these two equations, 4-9 and 4-10, are both of the form AX - B. Thus, 
this appears to be a hnear problem in X, but the emitter location solved by TDOA 
techniques is not linear. In this case, the fourth component of X is non-linear; this non- 
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linearity will act as a separate constraint on the solution of the linear problem. In the 
overdetermined case when there are four equations (M=5 observers), the non-linear 
constraint is automatically satisfied; however, we rarely expect to be in an overdetermined 
case. Therefore, the matrix A is singular, but still can be inverted as a generalized 
inverse; A~^ = A^(^AA^^ . The best way to obtain A'^, however, is by the singular 

value decomposition (SVD) algorithm. More compelling, the SVD algorithm directly 
provides the null space, , required to satisfy the non-linear constraint. Since A is 3x4, 
and we assume rank 3, then the null space of A is one dimensional and the complete 
solution for X is; 


X = A ^B + Aa^ 4-1 1 

Where is a unit vector spanning the null space of A and /I is a real variable. To 
prove this is the solution to AX = B , multiply both sides of the solution by A and recall 
Aa- = 0 , 


AX = AA^(AA^y'B + Xia^ = AX = B 4-12 

Given the above equation, it is evident that the constant A is arbitrary with respect 
to the linear equation; however, from equation 4-9, the non-linear constraint requires the 
fourth component of X be equal to the magnitude of the vector of its first three 
components; ||F'f - (A' 4 )' = 0 . 

Now, consider the first three dimensions or components of the solution for X — 
Chaffee and Abel call this vector; = A'{aA') 'u = A^^B. The solution to the linear 

equation has given us this vector which describes the emitter’s location uniquely in a 
three-dimensional subspace orthogonal to , but cannot provide any information to the 
ambiguous component normal to this subspace, i.e. along a" . The non-linear constraint 
must be used to determine this component defined by x.. Substituting the four 
dimensional solution to the linear equation, Chaffee and Abel label , plus Aa" back 

into the Lorentz functional, yields a quadratic equation forx : 

{a\a^)A^ +2A{w^,a^) + iw^,wy) = 0 4-13 
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which provides a direct solution for A. 

Thus, to summarize the result and combining the two notations; 


where 

= 


+ Aci^ ) — ^ 
~r'^ 


r 


n 


= w,+ Aa^ + Z 

I 

b 

A 

ji_ 


which is from the constraint 


4-14 

on the solution vector 



The quantity b in Equation 4-14 is called the bias which, used in GPS 


equations, represents the offset of time between the master’s and emitter’s time reference. 
The vector r gives the emitter’s location. 

Like all direct solutions, the impact of squaring the TO A equation results in two 
possible solutions for the emitter location. For GPS, this is not a problem since one 
solution is on the face of the earth and the other lies opposite to the satellite plane in 
space. In a tactical environment, however, this problem must be solved with some prior 
knowledge of the target of interest or with some common sense with respect to the 
components which describe the target’s location. 


C. LINEARIZED LEAST SQUARES APPROACH 

The linearized least squares approach to the geolocation of an emitter provides a 
statistically robust estimate of the emitter’s state and covariance uncertainty estimate. 
While the method does require an initial guess with respect to the location of the emitter, 
in many cases this can be provided by picking a point in the area of observers or by first 
solving for the state via one of the direct methods mentioned above. Certanly, the 
performance of linear least squares algorithms is well-known and stable givai that care is 
taken into their formulation. In this regard, JBSTAFSim utilizes a Square Root 
Information Filter (SRIF), see Appendix A for a detailed description, which has been 
proven to be more numerically stable and accurate than other methods. In addition, the 
SRIF allows for elegantly simple inclusion of prior estimates to be included with the 
current state estimate for batch processing as an extended Kalman Filter. Further a hybrid 
linearized least squares solution aided by a direct solution will most likely be the 
geolocation method of choice by the ARL:UT system. In the mixed TDOA/FDOA case. 
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JBSTAFSim also enhances the numerical accuracy of the SRIF solution by scaling the 
partials matrix and TDOA and FDOA data so that every equation is dimensionless and on 
unity order. Appendix B contains a detailed description of the scaling of mixed 
TDOA/FDOA data in order to improve the numerical accuracy and stability of a linearized 
least squares solution. 

Consider in detail the mixed TDOA/FDOA problem as described in Reference 19; 


Observation equations: 


c TDOA,j{r) = R^'^ 

Av^'^ 


AvO) 


4-15 

4-16 


Measurement covariance matrix: 


M = 


^ ^ TDOA 


It- 


c a 


FDOA. 


4-17 


Linearized least squares equation: 




0 



-To 


c TDOAy ^ 
— FDOA, 

U V 




4-18 


Direction cosines are denoted by a, and R is the observer-target separation. In 
the state vector, dx and dv refer to the offset from the position and velocity, respectively, 
of the nominal state, and the subscripts p and L refer to components parallel and 

perpendicular, respectively, to the line of sight. Thus, for mixed T/FDOA, the above 
partials matrix is 2(M - l)xN, where M equals the number of observers and N is the size 
of the unknown emitter’s state. When only TDOA information is used for geolocation, 
the partials matrix is reduced to (M - l)xN, where N < 3, and, thus, no information with 
respect to the target’s velocity can be given. For the purposes of this thesis, only cases 
which are evenly determined are considered. The overdetermined solutions will most 
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likely always dominate the evenly determined ones; further, this method of analysis 
assumes that resources are typically scarce and, thus, the overdetermined case event is 
unlikely to occur. With this philosophy in mind, JBSTAFSim is designed to produce two- 
or three-dimensional fixes in the TDOA only case and four- or six-dimensional fixes in the 
mixed T/FDO A mode. 

In addition to the solution to the state of the target, the linearized least squares 
approach also provides an estimate of the error covariance of the state estimate. Viewing 
4-18 as the set of equations Ax = b, the estimate of the state vector takes on the form 
^ = P. If we are interested in only the first n components of the solution 

vector and take the first nm components of P correspondingly; then, given the implicit 
assumptions of Normality in the least squares and extended Kalman filter, it can be shown 
that; 


- pf Pn' (^n ~ ~ Xl 4-19 

where // represents the true location of the target. Since this covariance matrix 

parameterized by n is always positive definite, it can be said that its column vectors form 
the basis of the eigenspace and the corresponding eigenvalues are: Thus a 

hyperellipsoid defined by; 




a, 


a-y 


yL 


=t 


4-20 


forms surfaces of equal probability densities. Therefore, for a certain £, the probability 
that a point will lie inside the ellipsoid can be computed. For example, if an n=2 
dimensional fix with an ^ = 3 sigma error ellipse is desired, the ellipse defined by the 
covariance matrix would yield an ellipse for which there would be a 99.7% probability that 
the target lies inside the ellipse. Thus, the size of the ellipsoid created by the covariance 
matrix can be thought of as a measure of the uncertainty of the solution. 

In addition to generating error ellipsoids, the covariance matrix also provides a 
means to measure the amount by which the observer geometry is degrading the fix. This 
measurement is often referred to as the geometric dilution of precision (GDOP). Again, 
given that the residuals of the state vector are Normally distributed, the covariance matrix 
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is a measure of their statistical properties. If these errors are small, a linear approximation 
which relates them to the observations may be considered satisfactory. “Moreover, if the 
observation residuals are Normally distributed, the fix residual residuals computed via the 
linear approximation are automatically Normally distributed.”[Ref 15] In this case the 
linear error propagation may be defined from Reference 15 as; 

C^ = AMA\ 4.21 

The total error on the fix as the root-mean-square deviation is 

^Tr[c^) = ^Tr(AMA^) 4-22 

which provides the “total fix position error.” Dividing this quantity, as suggested above, 
by the total measurement error provides an estimate of the relative dilution of precision: 

^Tr(AMA^) 4-23 

ylTr{M) 

While there is nothing geometric about 4-23, it places all the geometric errors in the 
{AA^) factor separately from the stochastic contribution of C^. The relationship 

between GDOP, measurement error, and fix error has been described as: 

^position GDOP ^measuremem 4-24 

Therefore, the GDOP is not truly geometric, but depends on the stochastic portion of the 
model. It is, however, “accurate for quantitative predictions within a class of stochastic 
models defined by covariance matrices that are not too different from each 
other”[Ref 15]. Further, it appears that based on the numerical experience of the 
ARL;UT test system, that it is a good tool for quantitative prediction.[Ref 20] As will be 
shown in Chapter VII, while GDOP provides a good relative measure of the geometry, a 
better means to measure the degree to which the geometry of the observers is needed. 
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V. JAVA-BASED STOCHASTIC TDOA AND FDOA SIMULATION 

(JBSTAFSIM) 

A. BACKGROUND 

JBSTAFSim was written to stochastically explore the differences and similarities 
of TDOA and TDOA/FDOA GPS-assisted geo-location solutions. It assumes a Multi¬ 
platform, GPS-assisted architecture similar to that being developed by ARL:UT as 
described in the previous chapters. JBSTAFSim allows the user to: create scenarios of 
various types of sensors and targets defined by straight-forward parameters; distribute 
these platforms in a simulated operating environment; dictate solver characteristics for the 
linear-least squares formulation; and adjust the measurement error distribution and 
correlation. While the purpose of JBSTAFSim was to mainly explore the differences 
between the two types of solutions in the evenly determined case, it could also be used a 
tool to investigate and experiment with various solver properties, solution methods, or 
algorithms in a variety of cases. Further, it could simply be used to experiment with 
various platform combinations, characteristics, or allocations to provide an understanding 
of the system’s strengths, weaknesses, and limitations. Finally, JBSTAFSim is an object- 
oriented program which has been written to be extendible. Given an understanding of the 
model as described in the following sections, the simulation can be extended to satisfy any 
of the areas of study discussed above or modified to include more complex sensor or 
target models. 

B. DISCRETE EVENT SIMULATION EN JBSTAFSIM 

If geo-location of an RF transmitter using TDOA and FDOA techniques was a 
straight-forward analytical problem, then a simple closed-form solution could be produced 
on paper. Yet, the signal processing techniques and assumptions presented in Chapter III, 
show that any realistic geo-location system is a “system” built on highly complex models 
so complex that any one analytic model designed to mathematically represent the whole 
system will fail to realistically describe the real-world system. Simulation, on the other 
hand, provides a means to “evaluate a model numerically, and ... [gather data] to estimate 
the desired true characteristics of the model” [Ref 23]. Therefore, simulation provides a 
means to examine the output measures of effectiveness by numerically exercising the 
model with inputs that would violate the assumptions of the analytic model and, thus. 
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provide estimates that are more characteristic of the real-world system. Given this 
complex mathematical model which acts dynamically and is influenced by stochastic 
components, JBSTAFSim, then, provides an evaluative, numeric tool which can evaluate 
and/or estimate the performance of a Muhi-platform GPS-assisted geo-location system as 
it would perform in a realistic environment—an analysis well outside and beyond the 
bounds of that provided by any mathematical formula, theorem, or set of equations. 

In Discrete Event Simulation (DES), simulated time only passes when events are 
not occurring. “In mathematical terms, we might say that the system can change [state] at 
only a countable number of points in time” [Ref 23] Events are listed on a master 
schedule in time and priority order. The entity which passes simulated time and ensures 
that events are executed in simulated time is often referred to as the Time Master. When 
no events are occurring, the Time Master simply advances the simulation clock to the 
scheduled time for the next event on the event list. The Time Master then allows that 
event to execute; this event may cause other events to execute or be scheduled with the 
Time Master. When the event is complete, the Time Master takes control again and 
executes the next event on the event list, advancing time as appropriate. When there are 
no more events left to schedule or execute, the simulation is finished. 

While Java is not strictly a simulation language, as Java gains popularity, new 
packages are being developed for many Java-based applications, including simulation. 
JBSTAFSim utilizes a simulation package named SEMKIT which was developed at Naval 
Postgraduate School by LT Kirk Stork, USN, in conjunction with Professor Arnold H. 
Buss. SIMKJT is a package designed to facilitate Discrete Event Simulation (DES). The 
version of SIMKIT used in JBSTAFSim is a development release written in the JDK 1.0.2 
and utilizes the version of Object Space’s Java Generic Library (JGL) written for the JDK 
1.0.2. Since SIMKIT is under development, further references to SIMKIT shall refer to 
the development release version of SIMKIT used in JBSTAFSim. Detailed information 
about the JGL and SIMKIT can be found at www.objectspace.com and 
http://131.120.142.115/~stork/simkit_home.html, respectively. 

In addition to providing the DES mechanization, SIMKIT also provides a 
“RandomStream” class to produce psuedorandom variables. While Java does implement a 
psuedorandom number generator, SIMKIT’s RandomStream class not only provides more 
suitably random numbers, but also provides the user a library of random variate 
distributions. JBSTAFSim employs several of these psuedorandom number generators for 
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platform movement and makes use of SIMKIT’s random uniform [0,1] generator to 
produce the random distributions for TDOA and FDOA measurement errors. 

C. MODEL ASSUMPTIONS 

Before going too deeply into model assumptions and sensor models, it is important 
to understand that JBSTAFSim was designed to provide a stochastic simulation for the 
comparison of TDOA and TDOA/FDOA geolocation systems. Like any good model, the 
model strives for simplicity wherever possible. As described below, platforms which move 
in the simulated environment are fairly simple yet retain their basic qualities which 
differentiate them from each other. Further, the logic of the solver is a very limited to 
solving for a target’s position given the data provided by the sensors. The solver expects 
an evenly determined case and it does not accept user guidance during the simulation nor 
does it make decisions using a complex algorithm or neural network. Functions such as 
these are those which are assumed to be in place in an operational system. In addition to 
assuming a geo-location system which is properly functioning and that the master has 
good communications with its slave sensors, the following sections describe the clarifying 
assumptions made in JBSTAFSim. 

1. GPS Data 

First, JBSTAFSim assumes that the platform GPS positions and velocities 
represent truth. For the model, this is a necessary condition since the model must establish 
some representation of ground truth. In reality, this is also a fairly good assumption given 
that a survey-grade GPS unit, which is fully Precise Positioning Service (PPS) and 
precision (p-code) capable, is implemented in the ARL;UT test platform. With this type of 
equipment, GPS fixes are considered to be accurate to within 10 meters. With respect to 
velocities, the manufacturer of the GPS receiver in the ARL;UT prototype system asserts 
velocity measurements accuracy of 0.01 meters per second. [Ref 24] 

As a consequence of using GPS information, the coordinate system used by 
JBSTAFSim and the ARL:UT system is the earth-centered-earth-fixed coordinate system 
as established by the World Geodetic System of 1984 (WGS84 ellipsoid). This system 
places all platforms in an ellipsoid defined by three coordinates (X, Y, Z) in units of 
meters. In order to assist the user in inputting data and for graphically displaying results. 
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JBSTAFSim receives and displays data in coordinates Latitude, Longitude, and Height. 
All calculations and results, however, are measured in meters in the WGS84 elUpsoid. 

2. Platform Models 

The platform models which function in the JBSTAFSim simulated environment 
require some clarifying assumptions for each platform’s movement within this simulated 
space. In order to understand how the JBSTAFSim platform models fimction in the DES 
environment. Figures 5-1 through 5-4 depict event graphs for each platform. 

CL Ship Model 

The ship model retains the most common features of all five platform 
models. Like all platform models, the ship keeps track of its current position and velocity 
in earth centered earth fixed coordinates. The event graph for a ship’s movement is shown 
in Figure 5-1. Upon creation, any ship in JBSTAFSim is provided by the user an 
operating box in which to steam. Also at creation, the ship chooses independent random 
uniform X and Y coordinates in this box as its destination. The ship calculates the amount 
of time it will take to arrive at its destination based on its current speed and distance to the 
destination. The ship then schedules an event to create a new destination when it arrives 
at the current destination. In addition to destination events, ships schedule their own 
change speed events. New ship speeds are selected randomly fi’om a uniform distribution 
defined by the low and high speed bounds provided by the user. Occurrences of speed 
change events are modeled as events with an exponential distribution with input parameter 
provided by the user. The first speed change event is scheduled upon initialization. When 
a speed change event occurs, in addition to changing the speed of the ship, the waiting 
change destination event is canceled and a revised arrival time is computed based on the 
new speed of the ship. A new change destination event is then scheduled for this revised 
time. Finally, the speed change event schedules a new speed change event to occur at a 
random time provided by the exponential distribution. 
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b. Ground Unit Model 

The ground unit model is slightly more complex than the ship model. The 
event graph for the ground unit model is provided in Figure 5-2. Like the ship model the 
ground unit keeps track of its position and velocity and schedules change destination 
events and change speed events. The difference with the ground unit is that it occasionally 
ceases to move, or does what is called a “take break event.” Break events are scheduled 
as a random process with an exponential distribution with input parameter provided by the 
user. The length of such take break events are also modeled randomly by an exponential 
distribution with input parameter provided by the user. When a take break event occurs, 
in addition to zeroing the speed of the ground unit, the waiting change speed event and 
change destination event are canceled A revised change destination event is scheduled for 
the remaining transit time to the current destination but only after the amount of break 
time has elapsed. Similarly a change speed event is scheduled in an amount of time 
provided by the random distribution plus the amount of rest time for the break. Finally, 
after the rest period has elapsed, a take break event schedules another take break event in 
a random amount of time provided from the exponential distribution. 
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Figure 5-2, Ground Unit Model Event Graph. 


c. Aircraft Model 


The aircraft model is more simple. Each aircraft patrols a circular “orbit” 
defined by the circle’s center and radius. In addition to these parameters, the user also 
defines the aircraft’s base height and how much the aircraft deviates from this base. Since 
the aircraft moves around its center in a fixed two-dimensional circle, the aircraft’s speed 
is vital to determining where on the circle the aircraft lies. For any circle with radius r the 
angle 6 at the center of the circle swept by an aircraft traveling along an arc length 5 is 

simply defined as 0 = —. The rate at which this angle changes can be simply defined by 


the first derivative: 


d6 1 ds ds 

— = —- = 0 ). Since — is a measure of the magnitude of the 
dt r dt dt * 


II I *1 

aircraft’s velocity, |v|, co can be computed as —. With this value of the aircraft’s 

position on the circle can be computed as a function of time. Thus, the aircraft model 
need only keep track of a and update its value when |v| changes. In order to prevent two 
aircraft from having the same or near similar values for 6 , each aircraft is given its own 
offset angle chosen randomly over the interval [0, 20], With this model, an aircraft need 


32 








only update its position on its circular track when polled and then pick a random height 
adjustment to make to its base height. However, when polled, the aircraft must also 
update its current velocity state. In order to compute the current velocity of an aircraft, 
the aircraft need only look a very short amount of time into the future and draw a straight 
line to this point to set its velocity vector. With this construct in mind, the event graph 
below is simple to understand. The aircraft need only update positions and velocity when 
polled or when a speed change event changes |v|. As before, speed change events occur 
randomly; are modeled with an exponential distribution; and cause a new speed change 
event to be scheduled. 



Figure 5-3. Aircraft Model Event Graph. 


d. Satellite Model 

While much could be done to model satellites of different orbital 
characteristics, the satellite model in JBSTAFSim is remarkably simple The reason for 
this simplicity is for the model to serve its purpose as an evaluative simulation. Certainly, 
more complex orbits could be included by extending the current model. JBSTAFSim’s 
satellite model assumes a satellite, or satellites, in a circular orbit with no gap in coverage. 
The model only requires the orbit height, starting and ending footprints of the region of 
coverage. Simulating continuous coverage, the satellite flies starting point to endpoint and 
returns back to the starting point. Velocity information need not be entered since, for 




satellites in a circular orbit, the magnitude of the satellite’s velocity is determined by the 

, , / GM „ w 

height of its orbit: |v| = J-—, where G is the Gravitational Constant in ^ Ag ^, M is the 

mass of the earth, r is the radius (height) of the orbit, and R is the radius of the earth. 
Since the satellite’s movement is deterministic, like the aircraft model, satellite position 
and velocity are updated only when polled. Otherwise, the satellite continues across the 
sky and returns to its starting point when it reaches the ending footprint point. 



Figure 5-4. Satellite Model Event Graph. 


e. Fixed Model 

The simplest of all models is the fixed sensor. The fixed model’s position 
never changes and its velocity is set to zero. A call to update a fixed object’s position 
literally results in nothing happening. Simply, the fixed model sends its information on 
whenever polled and requires no discrete event mechanization. 

3. Solver 

By using a square root information filter (SRIF) in a linearized least-squares 
formulation, JBSTAFSim’s solver assumes that measurement errors are statistically 
independent and distributed Normal with mean zero and input covariance. With some 
modification, a weighting matrix that assumes some correlation could be implemented; 
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however, since little is known about the true distribution of TDOA and FDOA 
measurement errors, an estimate of the correlation seems even more unlikely. 

In order to have an initial estimate of the target’s position, the linearized least 
squares geo-location solution is started with the results provided by a direct solution. 
JBSTAFSim solves this set of equations as described in Chapter IV using a Singular Value 
Decomposition routine adapted from Reference 21. As discussed in Chapter IV, after a 
fix is computed, the solver then uses the previous solution as the starting position for the 
subsequent solutions until a request is made for another direct solution to update the 
current solution. The periodicity of these updates may be set by the user. In an actual 
Geo-location System, direct solution updates might be requested as a function of many 
complex variables and parameters, for instance; quality of the previous fix, the 
dimensionless shock, time since the last fix, by request of the system user, or quality of the 
received signal. While having a set doctrine for computing direct solutions leaves little 
room for qualitative decision making algorithms that would exist in the ARL'.UT system, it 
is a simple abstraction for the simulation given the assumptions already made. In addition, 
by updating the target’s initial “guess” position, the direct solution prevents the model 
from losing track of the target if the previous linearized least squares solution was too 
great in error. The difficulty with any direct solution lies in the dual solutions provided as 
a result of squaring the TOA equations. As described earlier in Chapter IV, this ambiguity 
is generally resolved with some knowledge of where the target is in relation to the 
observers or by simply comparing the two solutions and throwing out the one that is 
obviously too distant from the observers. JBSTAFSim assumes that the operational 
system will have several algorithms for resolving this ambiguity, and, thus, always picks 
the correct starting point. 


4. Measurement Errors 

Related to the solver error assumptions, is the actual error mechanization itself 
JBSTAFSim measures precise TOA’s and FOA’s in meters and meters per second These 
values are differenced from those of the platform designated as the master. Thus, these 
differences are the true TDOA’s and FDOA’s that any sensor would measure in the 
absence of measurement noise. In order to simulate measurement noise, the truth values 
are then corrupted by Normal, zero mean random values generated from a distribution 
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defined by the user. The user is free to define the standard deviation of errors across each 
TDOA/FDOA measurement. In addition, the user is allowed to correlate the errors 
between TDOA/FDOA measurements. 

The production of measurement errors begins by first producing a covariance 
matrix using the data provided by the user. While the user is only concerned with the 
correlation of measurement errors between those differences which include the master as a 
sensor and not the correlation between slaves, for the three dimensional case, this 
correlation must nevertheless be supplied in order to build a covariance matrix from which 
to generate random errors. Further, JBSTAFSim checks this value against the partial 
correlation with those errors associated with the master to ensure that the given 
correlation vector generates a positive-definite symmetric covariance matrix from which to 
generate errors. The covariance matrix is then factored into a lower triangular matrix 
using the Cholesky decomposition routine. Next, JBSTAFSim uses SIMKIT’s uniform 
[0,1] generator to produce Normal (0,1) variables using the Box and Muller method. 
Calling the Normal (0,1) variables Z,, the correlated errors X^, and the lower triangular 

matrix C, the measurement errors are produced from the following algorithm (Scheuer 
and Stoller’s method) from Reference 23: 


For/= X,='^C,jZj 


./=! 


5-1 


Equation 5-1 could be changed to produce non-zero mean random numbers by 
simply adding the mean, //,, to the sum on the right hand side. Equation 5-1 also makes 

clear the need for a positive definite covariance matrix since any other matrix could not be 
factored into the matrix C. 

D. MODEL DATA PROCESSING 

As discussed above, JBSTAFSim assumes a GPS-assisted geo-location system that 
is functioning nominally with good communications between all sensors. JBSTAFSim 
sensors are polled and fixes are computed every 10 seconds of simulation time. While this 
poll frequency is not set permanently, it was chosen to correspond in value to that in some 
initial tests of the ARL:UT system JBSTAFSim never misses a fix or loses 
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communication between sensors since these are not characteristics that JBSTAFSim is 
endeavoring to simulate. 

The polling, fix, and result comparison process are depicted in the flow-chart in 
Figure 5-5 below. 



In addition to scheduling the next poll sensors event, the poll sensors event causes 
the target and all sensors to update their positions and velocities. Sensors, who maintain a 
reference to the target, compute their TO A (distance to target measured in meters) and 
FOA (line of sight velocity measured in meters per second) to the target. The slave 
measurements are then differenced from the master to form N-1 TDOA’s/FDOA’s. This 
“truth” data is then corrupted by inducing measurement errors that are modeled by a 
Normal distribution with zero mean and standard deviation and cross-correlation input by 
the user. This data and the sensor positions are then passed to the solver. The solver uses 
the sensor positions and the current guess of the target’s position to form the partials 
matrix for the square root information filter (SRIF). The SRJF then performs the required 
iterations to solve for the target’s position. These results are sent to a data collector that 
compares the solution to the truth position The data collector keeps running statistics— 
mean, variance, 95% confidence interval half-width, largest and smallest extremes—of the 
miss distance, geometric dilution of precision (GDOP), and size of the area of uncertainty. 
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Finally, the solver need not be a SRIF. Due to JBSTAFSim’s object oriented 
design, the data returned by all sensors to the master is itself an object which could be 
processed by any solver. This data object consists of a collection of objects which hold 
each sensor’s position, velocity, TOA data, FOA data, and has methods to evaluate the 
numerical quantities which make up the partials matrix. The master’s data is stored in the 
first location of this data object vector. In order to form the “DOA’s” from the data, the 
remaining data objects difference their information from the master’s. Of course, in a 
linearized least squares formulation, the partials matrix must be evaluated at some point; 
this point is supplied by an object which contains the current target estimate position and 
velocity. When the solver finishes its iterations, it returns its solution as this type of an 
object. Thus, any solver coupled to JBSTAFSim need only accept the basic objects in 
JBSTAFSim which define position and velocity for any platform. 

E. STATISTICAL RESULT ESTIMATION AND CALCULTAIONS 

Since the focus of analysis is centered on the miss distance of the geo-location fix, 
JBSTAFSim uses this parameter’s precision for the terminating condition. Miss distance 
is simply the Euclidean distance between two points—the truth position and the solution. 
Since both of these points are represented as objects in Java, the calculation of the miss 
distance is merely a matter of asking one object the distance to the other. 

Thus, JBSTAFSim seeks a precise estimate for the mean of the miss distance. 
However, since JBSTAFSim runs for no fixed sample size, the model has no control over 
the confidence-interval half-length which directly impacts the precision of the estimate. 
Therefore, in order to produce a precise estimate for the mean miss distance and provide a 
terminating condition, JBSTAFSim employs relative precision. Using relative precision as 
a terminating condition for the estimation of the mean of a random variable is discussed in 
Appendix C and in detail in Reference 23. The advantage of using relative precision lies in 
its ability to provide a precise estimate of the mean miss distance when the sample size is 
beyond the analyst’s control. 

In addition to fix miss distance, JBSTAFSim keeps statistical information on the 
geometric dilution of precision (GDOP) and uncertainty ellipsoid area or volume. 
Ellipsoid volume is calculated using the eigenspace as described in Chapter IV, Equations 
4-19 and 4-20. JBSTAFSim employs a Jacobi Transformation algorithm which consists of 
a sequence of orthogonal similarity transformations. Each transformation is a plane 
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rotation designed to annihilate one of the off-diagonal matrix elements. While each 
transformation may undo previously zeroed elements, each element continues to get 
smaller and smaller until the matrix is diagonal to machine precision. [Ref 21] With the 
eigenvalues of the estimate covariance matrix, the area or volume of the ellipsoid of 
uncertainty can be calculated as described in Chapter IV. GDOP information is 
accumulated to rate the quality and consistency of the geometry of the simulation. GDOP 
calculation is a straightforward application of Equation 4-23. 

These statistical results and accumulated data are provided to the user through the 
model output selections provided by the user at the simulation’s initialization. 

F. MODEL INPUT AND OUTPUT 


1. Input 

For any simulation run, JBSTAFSim requires a large amount of initialization data. 
The solver needs to know how many dimensions and how many sensors with which it will 
be dealing. In addition the user must set the weights for the TDOA and FDOA data. In 
practice, these weights would be computed on-line by the system using a hybrid of both 
equations and algorithms similar to those in Chapter 3 and from experience of the system’s 
performance. The solver can also be set to “batch” mode in which it combines the 
information of the previous fixes with the current data to produce a fix. In addition to 
solver information, the user must set the initialization data for each sensor and target—up 
to five total platforms. As discussed above, by inputting the error standard deviation and 
error cross-correlation between sensor measurements, the user must also create a 
distribution for the measurement errors. 

While this amount of input is extensive, once it is entered, the model takes over 
and requires no further guidance by the user. In order to make the model initialization 
process easier, JBSTAFSim allows the user to save the input data to a file specified by the 
user. With this file of initialization data, the JBSTAFSim simulation run may be 
immediately started by simply selecting the input file. With this construct, any 
initialization need only be entered once; further the data stored in the input file is relatively 
easy to edit and could be changed in any text editor. 
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2. Output 


As the user chooses, JBSTAFSim outputs simulation results to the screen and/or a 
data file. If output to a file is selected, JBSTAFSim returns the results of its simulation to 
a file designated by the user. This file contains the statistical results of three data 
accumulators for each solution type. As described above, the statistical accumulators 
keep track of the fix miss distance, GDOP, and size of the uncertainty ellipsoid. At 
simulation termination, the number of data points, mean, variance, standard deviation, 
largest and smallest extremes of the data for each accumulator are written to the output 
file. 

If a graphical display has been requested by the user, JBSTAFSim animates the 
simulated movement of sensors and the target while displaying the TDOA and 
TDOA/FDOA fixes. In addition a small summary at the bottom of the animation provides 
the current fix miss distance, an updated average miss distance, and GDOP of the current 
fix. As previously mentioned, the graphical display is in the latitude, longitude, height 
coordinate system—although height is not displayed in the animation. The latitude scale 
on the left of the screen provides a quick reference for the scale of the display. Each tenth 
of a degree of latitude represents six nautical miles or 12,000 yards. Thus, given that the 
coordinates of each sensor in the simulation is represented by WGS84 Ellipsoid 
coordinates on the order of 10^ - lO’ meters which translates to order lO' latitude and 
longitude degrees, and is then translated into integer screen pixels displayed on a small 
scale (large area) chart, it should come as no surprise that most platforms in the display do 
not tend to move with great speed, especially when fixes are being taken every 10 seconds 
of simulated time. 

Despite these facts and the fact that the screen images provide little analytical 
information, the graphical display does provide a very effective and quick means to 
ascertain what is actually happening in the model for both the analyst and for any audience 
of the analysis. Further, the ability to provide such an animated stochastic simulation 
across a variety of computer platforms must be considered a true asset for analyst who 
must present the analysis to an (any) audience. Yet, Java is able to deliver not only on this 
ability, but also allows the analyst and audience to be separated by thousands of miles with 
only the Internet and a Web browser between the two. Accordingly, the ability for 
JBSTAFSim to function as an applet is described in the next chapter. 
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VI. JBSTAFSBVI AND THE WORLD WIDE WEB 


As its name implies, JBSTAFSim is written in Java—specifically, the first release 
of the Java Development Kit (JDK 1.0.2). Java is an object-oriented programming 
language which is currently gaining increased popularity in the programming world. One 
of the reasons Java has become so popular is its run time library which gives Java 
bytecode platform independence. The same code (without recompiling) can be used on 
Windows 95, Solaris, UNIX, Macintosh, and so on. Given the exponential growth of the 
Internet, this is a true advantage, if not necessity, for any programming language. As a 
result of this portability, JBSTAFSim can not only run on any of the machines listed above 
as an application, but also can run as an applet on any machine using a Java enabled Web 
browser. The JBSTAFSim applet requires no special programming, code, recompilation, 
or platform. It only requires a Web browser, nothing else. Therefore, due to Java’s 
portability, JBSTAFSim can be run and its results examined instantly from anywhere and 
on any machine with little overhead. Thus, Java is a full featured programming language 
that provides the ability to be run from the Internet on any machine via a Web browser 
There are many more advantages to programming in Java, but simply stated, it is a simple, 
robust, portable, high performance, and dynamic language which allows for efficient, 
object-oriented, distributed code. For a detailed explanation of Java’s design and 
accomplishments. Sun has published a “White Paper” which can be found at 
‘‘httpy'/java.sun.com/whitePaper/java'-whUepaper-l.html. ” 

Java programs that work within a Web browser are called applets. The idea 
behind an applet is simple; “users download bytecodes from the Internet and run them on 
their own machines” [Ref 22]. Some might ask why use a Java applet on a Web page to 
run programs? Since Java is a full programming language, it has all the power of any 
“true” programming language with the ability to run dynamically on the host computer via 
the Internet. Without Java, if you were to design a “dynamic” Web page that could 
respond to a user or a stochastic simulation, each change in the Web page would require 
that data be sent to a CGI script on the server. The script would need to process the data 
and send the results back to the browser—probably requiring the creation of a new Web 
page on the fly. This is a horror to program and without a doubt will be slow. With Java, 
once the applet’s bytecodes are downloaded, the power of the programming language is 
available to run without updates from the server. In order to protect the host computer 
system from malicious bytecodes, Java implements a “sandbox” that restricts the applet 
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from corrupting memoiy outside the applet’s process space. In addition, secure Web 
browsers, like Netscape Navigator, or browsers that have their security features enabled 
by the user (Microsoft Internet Explorer) will prevent applets from writing or reading 
local files. Applets can include graphical user interfaces, can catch mouse movements and 
clicks, and can interpret text inputs. Further, all processing is performed by the user’s 
system; thus, the originating server is not continually bombarded with hits for information 
and number crunching. In addition, the user is not hampered by bandwidth when running 
an applet (although, bandwidth could be an issue dcwnloading the applet).[Ref 22] 

The power of Java and the applet concept above can be realized within the nature 
and abilities of the JBSTAFSim applet itself Quite simply, the only difference between the 
JBSTAFSim applet and JBSTAFSim application is one Java class which declares that it is 
an applet and will oversee the events which occur in the JBSTAFSim simulation. With no 
revisions or recompilations, every bit of code that the simulation uses as an application is 
used by the applet. Thus, the full power and stochastic nature of the JBSTAFSim 
application can be delivered anywhere to a Web browser via the Internet with precisely the 
same code. An illustration of the JBSTAFSim applet running in Netscape Navigator is 
shown in Figure 6-1. 

As an example of the power, flexibility and utility of this design, consider how the 
JBSTAFSim applet is currently configured to run from the World Wide Web. A user, 
somewhere on some operating system with a Java enabled Web browser starts the 
JBSTAFSim applet by visiting the appropriate Web page. This Web page currently exists 
on a server which is running on a SPARC clone with NEXT as its operating system. The 
JBSTAFSim bytecodes are retrieved by the server from another SPARC clone which uses 
SOLARIS as its operating system. The JBSTAFSim bytecodes consist of bytecodes 
compiled on a Windows-95 system. This code implements SIMKIT bytecodes which 
were compiled on a Macintosh. Further, the SIMKIT bytecodes implement bytecodes 
from the Java Generic Library that were compiled on some system unknown to the author 
of JBSTAFSim The server, ignorant to all of this, delivers all these bytecodes to the 
user’s system for execution in the user’s Web browser. The point is: it is irrelevant from 
what system the code originated; from where the code is coming; and to where the code is 
going for execution. The bytecodes need only be delivered to the user’s machine which 
will locally process the bytecodes’ information quickly and efficiently—all bytecodes run 
the same everywhere. 
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Figure 6-1 JBSTAFSim Applet in Netscape Navigator 


For both the analyst and the audience, this is a powerful tool. The analyst’s 
sponsors can quickly inspect a product or proposal without special software, hardware, or 
large investments of time and money. Further, the tool is extendible and can be coupled 
with other models regardless of the operating system or environment in which they were 
originally created. In the case of JBSTAFSim, in addition to ARL;UT and faculty at 
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Naval Postgraduate School, the applet’s audience quickly included the Naval Information 
Warfare Activity and other Commander, Naval Security Group offices and activities. 
Each activity only requiring a Web browser and a connection to the Internet. Certainly, in 
this light, JBSTAFSim could represent one small component in what has the potential to 
be a robust, high-performance, and distributed modeling architecture for the DOD and its 
community. 
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Vn. SIMULATION ANALYSIS 


A. BACKGROUND 

Several simulation scenarios were conducted to evaluate and compare the 
performance of both the TDOA and mixed TDOA/FDOA solutions. Both two 
dimensional and three dimensional fix scenarios were conducted. Within these two cases, 
several parameters were varied to quantify their effect on both solution types, including; 
measurement error variance, the correlation between measurement errors, geometry of 
the sensor network, and the frequency that the SRIF requests direct solution updates. 

B. TWO DIMENSIONAL GEO-LOCATION SCENARIO RESULTS 

The two dimensional case consisted of a sensor network of two ships and a 
satellite geo-locating a ship, as depicted in Figure 7-1. This scenario was run seven times 
each for four measurement error correlation values. The seven runs varied: the 
frequency—continuous, every 10 fixes, and every 20 fixes—of direct solution updates; the 
TDOA measurement error variance; and the FDOA measurement error variance. The 
measurement correlation values were: 0.0, 0.1, 0.5, and 0.9 error correlation between all 
measurements. The TDOA measurement error variance was examined for the values of 
60^ and 30^ meters. The FDOA measurement error variance was examined for the values 
of 0.2^ meters per second and 0.1^ meters per second. Sensor to target geometry was 
mediocre with GDOP averaging 1.5 with minimum and maximum values of 1.3 and 2.6, 
respectively. For reference, the values of 60^ meters and 0.2^ meters per second for the 
respective TDOA and FDOA measurement error variances are approximately equivalent 
to those empirically experienced in the ARL.UT system. 

1. Geo-Location Accuracy 

Results from these simulation runs are depicted in Figure 7-2 through Figure 7-17 
with some summarized information provided Table 7-1. The distance from the solution’s 
position to the target’s actual position—^the miss distance—is the primary measure of 
effectiveness (MOE), while the area of the 2-sigma uncertainty ellipse is a secondary 
MOE. Thus, the primary MOE is a measure of the solution’s accuracy while the 
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TD0A/7D0A Solucion 
fix Hiss Distance: 
41.8269 

Avg. Hiss Di3C4:nce; 

347.648 

GDOP: 

1.59271 


TDOA Solucion 
Fix Hiss Distance: 
41.8269 

Avg. Hiss Distance: 

472.293 

&DOP: 

1,59271 


Figure 7-1 Two Dimensional Scenario Sensor Allocation 


secondary MOE is a measure of the solution’s ability to rate its precision. In terms of the 
miss distance, the primary MOE, Figures 7-2 through 7-8 show that the mixed 
TDOA/FDOA solution maintains first order stochastic dominance for each simulation run 
which varied the fix update fi'equency in addition to the cases where the TDOA and 
FDOA measurement error variances were contrasted. 

a. Fix Updates 

In order to prevent JBSTAFSim’s solution fi'om diverging due to a poor 
target estimate, direct solutions are provided at a regular update interval. This interval 
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was varied from continuous updates, to updates after every ten fixes, and for updates after 
every 20 fixes. Clearly, as seen in Figure 7-2, the continuous direct solution updates are 
degrading the accuracy of both solutions, but particularly so for the TDOA solution. 
These continuous update distributions are so poor, they are dismissed from further 
analysis. Notice that while the TDOA solution and mixed T/FDOA solution are more 
similar in the “10 fix update” and “20 fix update” cases, the “10 fix update” scenario 
offers the best shape to this distribution. This can be seen in the relative location of the 
point in the distribution where the curve makes its “right hand turn” in the upper quantiles. 
Notice that this location tends to be higher and more to the left in the “10 fix update” 


scenano. 


Comparison of TDOA and FDOA Miss Distance Quantiles 
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Figure 7-2 Comparison of TDOA and T/FDOA Miss Distance Quantiles 
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Figure 7-3 Comparison of TDOA and T/FDOA Miss Distance Quantiles 


Comparison of TDOA and FDOA Miss Distance Quantiles 



Miss Distance Quantiles in Meters 

_ 10 Fix Updates. 0.1 Coir., Sigma TDOA = 60m. Sigma FDOA = 0 2 m/s _ 

Figure 7-4 Comparison of TDOA and T/FDOA Miss Distance Quantiles 
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Figure 7-5 Comparison of TDOA and T/FDOA Miss Distance Quantiles 


Comparison of TDOA and FDOA Miss Distance Quantiles 



Miss Distance Quantiles in Meters 

_ 20 Fix Updates, 0.1 Corr.. Sigma TDOA = 60m, Sigma FDOA = 0.1 m/s _ 

Figure 7-6 Comparison of TDOA and T/FDOA Miss Distance Quantiles 
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Figure 7-7 Comparison of TDOA and T/FDOA Miss Distance Quantiles 


Comparison of TDOA and FDOA Miss Distance Quantiles 



Miss Distance Quantiles in Meters 

_ 20 Fix Updates, 0 5 Corr , Sigma TDOA = 3Qm, Sigma FDOA = 02 m/s _ 

Figure 7-8 Comparison of TDOA and T/FDOA Miss Distance Quantiles 


50 




























b. Solution Type Peiformance 

In each scenario above, the mixed T/FDOA case dominates the TDOA 
solution. The T/FDOA solution’s dominance is largest in the upper portion of the 
distribution after the characteristic curve as noted above. In order to show these contrasts 
more clearly. Figure 7-9, below, plots some of the scenario results together In this figure, 
a horizontal line representing the 0.95^*' quantile is drawn to intersect the solution 
distribution curves This line intersects the mixed T/FDOA solution, with measurement 
error variances of 30^ m and 0.2^ m/s, 100 meters to the left of the TDOA solution with a 
measurement error variance of 30^ m. Thus, the is a significant difference between the 
two solution types in the upper areas of their distributions. 


Comparison of TDOA and T/FDOA Miss Distance Quantiles 



Miss Distance Quantiles in Meters 

_ 10 Fix Updates Except * _ 

Figure 7-9 Comparison of Several TDOA and T/FDOA Miss Distance Quantiles 

Even more interesting, notice that the mixed solution curve identified 
above (error variances of 30^ m and 0.2^ m/s) outperforms a mixed solution with an 
improved FDOA error variance of 0.1^ m/s and a degraded TDOA measurement error 
variance of 60^ m. Further, note that at this quantile, this second mixed solution cur\'e 
(error variances of 60" m and 0.1^ m/s) has started to dominate the TDOA solution \\dth 
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an improved TDOA measurement error variance of 30^ m. Thus, this T/FDOA solution 
outperforms that of the TDOA solution in the upper quantiles where either measurement 
error or geometry is degrading the accuracy of the fix. Therefore, despite a larger error 
variance in the TDOA measurement errors, a mixed T/FDOA solution provides more 
accurate fixes when facing situations that are more difficult for the TDOA solution to 
solve by itself Thus, the designers or users of a system such as this would desire a larger 
reduction the TDOA measurement error variance as compared to that in the FDOA 
measurement even though it is the FDOA information which will improve the accuracy of 
the fix. 


c. Upper Quantile Performance 

The most striking feature of these curves is the sharp, near right-angle turn 
that the distributions take at the upper quantiles. The effect of rotating any of the quantile 
plots 90 degrees to the left so that the actual quantiles are on the X-axis make the “right- 
angle turn” feature of the distribution appear even more pronounced. With this 
understanding of the shape of the distribution, gaining an intuition of the difference 
between the two solutions that transcends first order stochastic dominance can be gleaned. 
As pointed out, the area in which the two solution distributions differ the most can be seen 
in the separation between the curves in the higher quantiles. In this light, the mixed 
T/FDOA solution’s dominance can be applied to the outliers of the distribution. In the 
cases where noise or geometry degrade the TDOA solution, the additional information 
from the FDOA measurements increases the accuracy of the T/FDOA solution over that 
of TDOA solution. Thus, it is especially in these cases of bad geometry or large noise that 
the mixed solution outperforms the TDOA only solution. Unfortunately, at each 
individual data point, it is unclear from these results whether it is the geometry or the 
measurement error which causes the size of the error in the location of the target. Two 
interesting plots of the mixed T/FDOA performance distribution are shown in Figures 
7-10 and 7-11 below. Notice in Figure 7-10 that the lower FDOA measurement error 
variance leads to better results as would be expected. In addition, holding the FDOA 
measurement error variance at 0.2‘ m/s and decreasing the TDOA measurement error 
variance leads to a distribution which dominates the other two distributions. These are not 
unexpected results; however. Figure 7-11 turns these results around. The distribution 
with the worst measurement error variance dominates Notice that Figure 7-10 displays 
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T/FDOA Miss Distance Quantiles by Sigma TDOA and FDOA 



Miss Distance Quantiles in Meters 

_ 10 Fix Updates, 0.9 Corr., Sigma TDOA - 6Qnn* _ 

Figure 7-10 Comparison of TDOA and T/FDOA Miss Distance Quantiles 


Comparison of T/FDOA Miss Distance Quantiles by Sigma FDOA 



Miss Distance Quantiles In Meters 

_ 10 Fix Updates, 0.0 Corr., Sigma TDOA = 60m* _ 

Figure 7-11 Comparison of TDOA and T/FDOA Miss Distance Quantiles 







































distributions where the measurement error correlation is high—0.9; yet. Figure 7-11 
portrays the results from 0.0 measurement error correlation. It is difficult to describe why 
measurements with higher error correlation should outperform those with no correlation. 
Further, it is not obvious that measurements with a larger error distribution should 
outperform those with better measurement data. Yet, this is exactly what these figures 
indicate. These results lead to two avenues of discussion—^geometry and uncertainty. 

If the results of Figure 7-10 and 7-11 were functions of geometry, it is not 
clear from the GDOP measurements the degree to which geometry is degrading the 
accuracy of the fix. Recall from Chapter IV that GDOP is not really a measure of 
geometry but a description of the diagonal elements of the estimate covariance and 
weighting matrices. It is more a measure related to what the system thinks might happen 
versus that of what is actually happening in the system. In Reference 13, Chaffee and 
Abel, the authors of JBSTAFSim’s direct solution theory, discuss how the null space of 
the direct solution relates to the two roots of the direct solution. Perhaps, the null space 
could provide some information with regard to the geometry of the observers and their 
target, however this type of analysis might only describe the geometry of the TDOA 
surfaces. It is more than conceivable that the ability of the mixed T/FDOA solution to be 
more accurate in the highest quantiles of the miss distance distribution—the outliers—can 
be attributed to cases where the geometry of the TDOA surfaces is poor yet the geometry 
of the FDOA surfaces is sufficient enough to improve the mixed solution. In view of the 
shape of the miss distance distributions in the upper quantiles, this theory could describe 
the sharp “right hand turn” as described above. The solution to rating the geometry of a 
network of observers in relation to their target is beyond the scope of this thesis; however, 
JBSTAFSim would be the perfect tool to empirically experiment with this theme of 
geometry. The second issue, fix uncertainty, deserves a section unto itself to describe. 

2. Uncertainty of Geo-locations 

If geometry is not all or only part of the answer to the unexplained improvement in 
the accuracy of the solution in the presence of correlation, then, perhaps, the uncertainty 
of the estimate may provide another avenue to help explore this issue. Figures 7-12 
through 7-17 display the distribution of the area of the two-sigma uncertainty ellipse 
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associated with both solution types for the continuous, 10 fix, and 20 fix update scenarios. 
In these graphs, the standard span (the area inside the “whiskers”) of the distribution is 
drawn at 1.5 times the Inter-Quartile Range beyond the quartiles. All points outside the 
standard span are drawn individually. With this construction, it is easy to view the density 
of the outliers. 

When the TDOA solution is compared with that of the mixed T/FDOA solution, 
the main difference between the two plots is the density of the outliers. Another way of 
indicating the larger spread of the T/FDOA uncertainty ellipse size is through the 
comparison of the standard deviation of this area. These quantities are summarized in 
Table 7-1 below. 


Error 

Continuous Update 

10 Fix Update 

20 Fix Update 

Correla¬ 

Sigma 

Sigma 

Sigma 

Sigma 

Sigma 

Sigma 

tion. 

TDOA 

T/FDOA 

TDOA 

T/FDOA 

TDOA 

T/FDOA 

0.0 

13571.4 

56457.7 

5618.3 

13197.1 

4711 

10962 

0.1 

13175.6 

55483.4 

5875.9 

12765.0 

4901 

9659 

0.5 

13389.5 

55616.6 

5643.1 

14718.5 

4703 

9658 

0.9 

13708.7 

55613.3 

5576.7 

14228.6 

4610 

11007 


Table 7-1 Distribution of TDOA and T/FDOA Ellipse Size Standard Deviations 


First, notice that the ellipse size standard deviations for the mixed T/FDOA 
solutions are much larger than those of the TDOA solutions in all cases. This is confirmed 
by the outlier spread and densities in the graphs below. Second, notice that in the 10 Fix 
Update scenario (where the accuracy of the 0.0 correlation was compared to that of the 
0.9 correlation in Figures 7-10 and 7-11 above) the solution which produced more 
accurate results—0.9 correlation—has a larger standard deviation of ellipse area. 
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Figure 7-12 Comparison of TDOA Ellipse Areas by Correlation 


TDOA/FDOA 2-Sigma Ellipse Area Distributions by Meas. Error Corr. 
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0.1 m/s 


Figure 7-13 Comparison of T/FDOA Ellipse Areas by Correlation 
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TDOA 2-Sigma Ellipse Area Distributions by Meas. Error Corr. 



Measurement Errror Correlation 

_ 10 Fix Upd. Sigma TDOA = 60m _ 

Figure 7-14 Comparison of TDOA Ellipse Areas by Correlation 


TDOA/FDOA 2-Sigma Ellipse Area Distributions by Meas. Error Corr. 
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Figure 7-15 Comparison of T/FDOA Ellipse Areas by Correlation 
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TDOA 2-Sigma Ellipse Area Distributions by Meas. Error Corr. 



Measurement Enror Correlation 

_ 20 Fix Upd.. Sigma TDOA = 6Qm _ 

Figure 7-16 Comparison of TDOA Ellipse Areas by Correlation 


TDOA/FDOA 2-Sigma Ellipse Area Distributions by Meas. Error Corr. 



Measurement Errror Correlation 

_ 20 Fix Upd., Sigma TDOA = 6Qm. Sigma FDOA = 0.1 m/s _ 

Figure 7-17 Comparison of T/FDOA Ellipse Areas by Correlation 
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Now it is possible to tie geometry, error correlation, accuracy, and uncertainty 
together. Clearly, the T/FDOA solution dominates in accuracy; however it tends to have a 
larger spread of uncertainty with regards to its results. Yet, comparing the median and 
inter-quartile ranges of the uncertainty ellipse sizes for the two solutions, they are very 
similar. Thus, the mixed T/FDOA solution is able to transcend errors induced by 
geometry and measurement errors, possibly induced by correlation, and, thus, provide 
more accurate fixes—especially in the upper quantiles of the accuracy distributions. The 
price for such accuracy is reflected in the uncertainty of the estimate. This should not be 
upsetting, and, in fact, proves the robust nature of the Square Root Information Filter. 
Since the SREF reports each solution with no idea with regard to its accuracy (accuracy 
requires knowledge of the true target location), it can only rate each fix by the solution 
estimate covariance. Therefore, even though the SRJF provides a more accurate fix, it 
correctly reports that measurement errors, error correlation, or geometry are degrading its 
confidence in the solution that it has provided. 

C. THREE DIMENSIONAL GEO-LOCATION SCENARIO RESULTS 

The three dimensional scenarios included a base line “good” geometry scenario, to 
show the effects of error correlation, and the comparison of a “good” and “bad” 
geometry cases to show the effects of sensor allocation on geometry. No dominance of 
either solution type for the three-dimensional case is shown 

Setting up scenarios to compute three dimensional fixes is very difficult. This 
difficulty is due mostly to the geometry of the scenario. A three dimensional fix requires 
four observers in order to form three differences for the evenly determined case. The 
addition of the height dimension, of course, expands both the partials matrix (from a 2x2 
to a 3x3 matrix for the TDOA solution and from a 4x4 to a 6x6 matrix for the mixed 
T/FDOA solution) and the data vector (from 2x1 to 3x1 and from 4x1 to 6x1 for TDOA 
and T/FDOA formulations, respectively) in the linearized least squares formulation. As a 
result each difference pair between slave and master must allocate its position information 
from the partials matrix to its respective TDOA and/or FDOA measurement(s). The 
difficulty in geometry arises from the formulation of partial matrices from these difference 
pairs that are near singular. While this is straightforward in theory, it is difficult to explain 
concisely to a warrior allocating his or her assets to perform a passive geo-location. 
Further, the possibility of the construction of a near singular matrix is not as readily 
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apparent when viewing the relative three-dimensional geometry as compared to that of the 
two-dimensional geometry. Finally, and most profoundly for an operational system, this 
problem is further exacerbated when one considers that in practice, the warrior may have 
little information with regard to the target’s true position and, thus, defining the geometry 
a priori becomes an even more difficult task. Yet, this does not mean three-dimensional 
fixes need be abandoned. The scenarios which depict “good geometry” provide excellent 
fixes—many with median miss distances well below 100 meters—that outperform poor 
geometry two-dimensional fixes. The main difficulty with a three-dimensional fix, then, 
may be providing the warrior with the means to optimally allocate his or her assets. 

1. Good Geometry Results 

The good geometry scenario consisted of a sensor network of two ships, an 
aircraft, and a fixed sensor targeting a ship. A visual depiction of the scenario is provided 
in figure 7-18. This scenario was run for 20 measurement error correlation values from 
-0.9 to 0.95. For the purposes of identifying the correlation combinations, a simulation 
run labeled with a negative correlation value actually means that within the three TDOA 
measurements and within the three FDOA measurements, two measurements are 
correlated negatively by a magnitude equal to the positive correlation of the other pair of 
measurements. For example, if measurements 1 and 2 were negatively correlated, then 
measurements 1 and 3 would be positively correlated while measurements 2 and 3 would 
be negatively correlated. 

The results of this scenario are multifaceted. First, restricting the sensors and the 
target into very tight boxes ensured that the geometry of the TDOA measurements would 
be almost always perfect. Thus, the TDOA measurements dominated the solution with 
regard to the position of the target. Therefore, the identical results of both the TDOA and 
T/FDOA solutions proves the stability and accuracy of the mixed solution’s square root 
information filter—a particularly gratifying result given the formulation, scaling, and 
un-scaling of the linearized least squares equation as described in Chapter IV and 
Appendix B 

The surprising result, as depicted in Figure 7-19, is the improvement in the 
accuracy of the solution as the correlation between measurement errors increases. 
However, upon inspecting Figure 7-20, the opposite is true of the distribution of the 
volume of the uncertainty ellipses. Given the results from the two-dimensional case 
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TDOAyiDOA Solution 
Fix Hiss Distfluice; 
113.983 

Avg. Hiss Distance! 

91.5161 

GDOP; 

0.954033 


TDOA Solution 
Fix Hiss Distance: 
113,988 

Avg. Hiss Distance; 

91.5161 

GDOP: 

0.954038 


Figure 7-18 Geometry of Basic Scenario 


discussed above, this is no longer a surprise. While the SRIF has no idea that it is 
producing more accurate results in the highly correlated noise environment, it is robust 
enough to realize that the measurement errors which it is encountering are deviating from 
what the weighting matrix had told it to expect. The SRIF produces the best fix it can; 
yet, it honestly reports, in the uncertainty of the fix, that the measurement errors are not 
indicative of what it had expected 
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TDOA and T/FDOA Miss Distance Distributions by Meas. Error Correlatioi 



Measurement Errror Correlation 

_ Sigma Tdoa = 60m, Sigma FDOA ~ 0.2nn/s _ 

Figure 7-19 Comparison of Miss Distance Distributions by Correlation 



Figure 7-20 Comparison of Ellipse Volumes by Correlation 
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Comparison of TDOA and T/FDOA Miss Distance Quantiles by Error Corr 



Miss Distance Quantiles in Meters 

_ Sigma TDOA = 60m. Sigma FPOA = 0.2 m/s _ 

Figure 7-21 Comparison of Miss Distance Distributions by Correlation 

2. Comparison of Good and Bad Geometry Scenarios 

As described in the introduction to this section, the geometry of the observers 
relative to their target is crucial for an accurate three-dimensional solution. Figures 7-22 
and 7-23 show the configurations of the “bad geometry” and “good geometry” scenarios 
for this illustration, respectively. Figures 7-24 and 7-25 show the quantile distributions of 
these two scenarios. 

Examining the distribution of sensor platforms in Figure 7-22, it is not clear—even 
with the position of the target known—that this allocation of sensor assets would lead to a 
“bad geometry” scenario; yet, the GDOP in this scenario averaged 222 with minimum and 
maximum values of 149 and 1472, respectively. Compare this to the “good geometry 
scenario in Figure 7-23 whose GDOP averaged 0.8962 with minimum and maximum 
values of 0 8874 and 0.91277, respectively. Without this data from the simulation runs, it 
would have been impossible to predict these results. Neither case appears exceptionally 
threatening from the geometry standpoint; yet, it will be the scenario’s geometry that 
directly impacts the accuracy of the geo-location solution. 


63 















































































































































Lattud 



Java-based Stochastic TOOA/FOOA Simulation 


33.4 


33.3 


32.9 


32.6 


- 119.6 * 119.4 - 119.2 - 119.0 - 118.8 -118 6 - 118.4 - 118.2 

Longitude 


V * 




TDOAyfDOA Soiucion 
fix Hiss Distance: 

16229.9 

Avg. Hiss Distance: 

14348.9 
GDOP: 

192.446 


TDOA Solution 
fix Hiss Distance: 

13942.9 

Iliss Distance: 
14311.3 
GDOP: 

194.88S 


Figure 7-22 Bad Geometry Case 
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Figure 7-23 Good Geometry Case 














Bad Geometry TDOA and T/FDOA Miss Distance Quantiles 



Miss Distance Quantiles in Meters 
Sigma TDOA — 60nn, Sigma FDOA - 0.2 nVs, 0.0 Corr. 

Figure 7-25 Bad Geometry Miss Distance Distribution 


Upon examining Figures 7-24 and 7-25, the possible impact of poor geometry on 
the accuracy of the three-dimensional solution becomes clear. An explanation of what 
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makes the “good geometry” case “good” is required. In terms of sensors, the difference 
between the cases is the replacement of one of the aircraft observers from the “bad” case 
with a satellite in the “good” case. Given its altitude, the satellite is able to create a better 
three-dimensional surface when its position is differenced from that of the master As a 
result, the satellite “sees” height better than that of the aircraft who may as well be on the 
surface of the earth in comparison to the satellite’s height. The partials matrix, then, 
contains more useful information with regard to the third dimension. Thus, the surface 
defined by the difference of the satellite’s and the master’s information is better suited to 
complement the other surfaces formed by the other observers and, therefore, can provide a 
more accurate fix with this “good geometry” of measurement surfaces. 

The second obvious feature when examining Figures 7-24 and 7-25, is that, like 
the first good geometry case in the previous section, both the TDOA and T/FDOA miss 
distance distributions are identical. In this analysis, the author was not able to produce a 
case where the T/FDOA solution stochastically dominated the TDOA solution. While the 
construction of such a case may be possible, performance of both solutions would most 
likely be highly dominated by the quality of sensor geometry versus that of sensor 
measurements. Further, while individual cases can be shown where a single TDOA or 
mixed T/FDOA solution was obviously better than that of the other, the author was not 
able to show that this lead to the stochastic dominance of either solution type. Thus, the 
solution to improving three-dimensional fixes does not appear to be in the augmentation of 
FDOA information; however, FDOA information did not on the whole appear to degrade 
the distribution of the three-dimensional solution either. The consequences from these 
results lead to two avenues of discussion: solution formulation and sensor selection. 

A solution to improving the three-dimensional fix may be simply collapsing to two- 
dimensions. This analysis has limited itself to the exploration of solutions in the evenly 
determined case. Yet, it is likely that most overdetermined two-dimensional solutions will 
dominate or at least do as well as those in the evenly determined case. In the 
overdetermined case, it has not been shown to what extent FDOA information could 
improve or degrade the solution. Likewise, the quality of a three-dimensional fix from an 
overdetermined sensor network with more than four observers has not been shown Yet, 
having a plethora of sensors to distribute at will seems overly optimistic and would not 
suit the warrior who desires to know a lower bound on the performance of his or her 
equipment. Yet, these are both important analysises which should be performed to 
understand the operational capabilities of the true system 
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On the other hand, if the utility for a three-dimensional fix is high, then the warrior 
who needs such high quality information should be concerned with the selection of robust 
three-dimensional assets. In particular, in this analysis it has been shown that the satellite 
is an excellent sensor to assist in observing height. Yet, it was also seen in two- 
dimensional simulation runs that in the rare case when a satellite passes directly overhead 
the target, it can say little about the target’s two dimensional location since the satellite 
only “sees” height. Again, the proper allocation of assets sometimes may be out of the 
warrior’s control. Given that satellite coverage may not always be available, another 
three-dimensionally robust-sensor might be one of the prototype unmanned aerial vehicles 
(UAV’s) that is designed to operate at very high altitudes for long periods of time—for 
example, the Tier 11+ Conventional High Altitude Endurance (HAE) UAV (Global Hawk) 
or Tier Ill-Low Observable HAE UAV (Dark Star) both of which advertise surveillance 
altitudes greater than 18 kilometers and 12 kilometers, respectively. Note that the first 
three-dimensional “good geometry” scenario only used one aircraft. Further, both UAV’s 
and satellites may be less obtrusive ways of improving the sensor to target geometry when 
airspace restrictions or the requirement for the concealment of sensors dictate that the 
traditional sensors would not suffice. 

Like most decisions a military commander must make, the selection and allocation 
of sensors to perform passive geo-location is not a straightforward choice. Yet, further 
analysises of like those suggested above could assist the warrior in this decision process 
Again, a good tool for exploring any of these issues would be a slightly modified version 
of JBSTAFSim. 
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vm. CONCLUSIONS AND RECOMMENDATIONS 


A. CONCLUSIONS 

A simulation based methodology was used to explore the differences between geo¬ 
location solutions computed using both time-diflference-of-arrival (TDOA) and mixed 
TDOA and frequency-difference-of-arrival (FDOA) techniques using a Multi-platform, 
GPS-assisted architecture. This methodology assumed no one method to measure either 
the TDOA or FDOA; yet, allowed for the stochastic modeling of the errors in such 
measurements in such a general way that they could apply to any TDOA or FDOA 
measurement process. Likewise, the platforms which measure such quantities were 
modeled in abstract ways, yet retained each platform’s most basic behavior. In addition, 
the simulation portrayed joint platforms operating together in a joint operating 
environment. To enhance the numerical veracity of the solutions, the simulation 
implemented a numerically stable Square Root Information Filter (SRIF) which received 
initial target estimates based on a numerically stable direct solution implementing the 
Singular Value Decomposition routine. These numerically robust methods were used to 
provide the most numerically stable and accurate solutions possible in the face of poor 
sensor geometry and computer round-oflf error. Data taken during the simulation 
included: the solution miss distance from truth, the area or volume of the 2-sigma ellipse 
centered around the solution, and the geometric dilution of precision (GDOP). The 
analysis focused on the accuracy of the solutions represented by the fix miss distance, and 
the solution uncertainty, represented by the volume or area of the 2-sigma ellipse around 
the fix. 


1. Two Dimensional Conclusions 

In the two-dimensional problem, it was shown that a mixed T/FDOA solution 
maintains first order stochastic dominance over the TDOA solution in mediocre geometry. 
The price the mixed solution pays for this increased accuracy is a large increase in the 
spread of its fix uncertainty; however, both solutions maintained roughly the same 
partitioning of fix uncertainty around the center of their respective uncertainty 
distributions. Therefore, increasing the spread of the uncertainty distribution remains a 
small price to pay for the increased accuracy of the solution. The largest difference 
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between fix accuracy was seen in the upper quantiles of miss distances where the mixed 
T/FDOA solution significantly improved the solution fix accuracy over that of the TDOA 
solution. Thus, the FDOA surfaces appear to substantially improve the solution accuracy, 
especially when the TDOA surfaces would lead to poor fixes. 

In addition to the contrasts between solution types, it was also shown that while 
the direct solution prevents the linearized least squares solution fi'om diverging, the 
amount to which the SRIF was updated by the direct solution significantly impacted the 
accuracy of both solution types. In particular, solutions which used the direct solution as 
an initial estimate for every fix were so poor, their results were disregarded. 

2. Three Dimensional Conclusions 

As an exploration into the impact of correlated errors, measurement error 
correlation was varied across 20 values in a good geometry scenario. While both 
solutions produced more accurate results in the higher, positively correlated cases, the 
solutions also indicated more uncertainty in these fixes. Again, this confirmed the robust 
nature of the SRIF formulation. These results are useful since the actual system would 
implement sensors taking measurements in a common operating environment with 
common noise characteristics, thus, most likely creating similar measurement errors. 

The three-dimensional simulations did not show stochastic dominance of either 
solution type. While some simulations showed differences in the accuracy of individual 
solution type fixes, neither solution type was shown to stochastically dominate the other. 
Results which produced identical accuracy distributions for both the TDOA and mixed 
T/FDOA solutions across several scenarios proved the stability of the algorithm which 
solves the mixed T/FDOA solution. The inference drawn fi'om these simulations was that 
geometry, more than solution and measurement type, impacts the accuracy of the three- 
dimensional fix solution. Therefore, the choice of sensors which can more optimally 
measure three-dimensional information should be the first concern for those who require 
such fixes. In particular, it was shown that sensors with high altitude characteristics are 
well-suited for this task. 
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B. RECOMMENDATIONS 


1. System Use and Design 

The Multi-Platform, GPS-Assisted TDOA Geo-location System as designed by the 
Applied Research Laboratories, University of Texas, Austin is a robust and potentially 
highly functional passive geo-location system for both the Navy and the Department of 
Defense. FDOA information can be safely added to the geo-location formulation to 
augment the TDOA information for the solution of the target’s location. In the two- 
dimensional solution case, the FDOA information can significantly improve the accuracy 
of the geo-location. In the three dimensional case, the additional information does not 
appear to degrade the system’s performance in any significant manner. 

The linearized least squares formulation in the SRIF appears to be a very stable 
and robust formulation for solving the geo-location problem. The system which 
implements this formulation should implement some direct solution or other means to 
provide an roughly accurate estimate of the target’s location to start the SRIF’s solution 
iterations. While the actual system would not be computing the large number of fixes for 
one target that were required in this stochastic analysis, a means to periodically update the 
SRIF with a new target estimate should be implemented. Further, the actual system would 
need to request this information based on signal characteristics, computation of the 
dimensionless shock, or some combination of both solver and signal properties. In 
addition, such an algorithm might be able to re-weight the partials matrix in the presence 
of correlated noise. The ability to chose a more optimal measurement covariance matrix 
would significantly improve the fix estimate covariance and, thus, decrease the uncertainty 
of the estimate of the target’s location. 

An avenue to improving three-dimensional fixes was shown to be through the use 
of platforms which can measure height information more reliably. This consequence leads 
to two recommendations. First, the ARL:UT system should be incorporated and tested on 
satellite systems and as a potential payload for prototype UAV platforms. Second, since 
this matter is at its crux a geometry issue, a new method to rate the geometry of the sensor 
network, other than GDOP, should be pursued. For the warrior, in particular, such a 
method should be able to assist in rating the sensor network’s geometry given possible 
locations of the target. Such a method would allow the commander to allocate his or her 
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assets in an pseudo-optimal way given either possible locations of the target or given no a 
priori information with regard to the target’s location. 

2. Java-Based Stochastic TDOA and FDOA Simulation 

JBSTAFSim proved itself as a robust stochastic simulation for the analyses 
presented in Chapter VII. In order to study any of the proposed areas above, modified 
versions of JBSTAFSim would provide the perfect tool for the analysis of any of these 
issues. In addition, JBSTAFSim could be modified so that specific platforms or platform 
characteristics, for example a particular type of satellite orbit, could be evaluated for its 
ability to contribute to the geo-location solution. In addition, due to time and scope 
limitations, JBSTAFSim’s “batch mode” option for its solver was never stochastically 
tested with a static target. JBSTAFSim requires more programming so that it can make 
qualitative decisions with regard to the fi'equency with which it resets its batch mode 
solver. 

Investigating solution types in both the overdetermined two-dimensional and three- 
dimensional cases should be undertaken. While it is likely that these solutions would 
dominate the evenly determined ones, it is not clear to what extent FDOA information 
might affect these solutions. This would be a particularly interesting study given that 
overdetermined, “bad geometry’ three-dimensional fixes could potentially be collapsed to 
more accurate overdetermined two-dimensional fixes. With some effort JBSTAFSim 
could be modified to analyze this issues. 

Next, JBSTAFSim has modeled both the TDOA and FDOA measurement errors 
as Normally distributed random variables. The Normal distribution was implemented for 
its ease in the ability to generate correlated random variables on a computer and, as 
presented in Chapter III, due to the fact that any measurement error distribution is in itself 
a fiinction of many random variables which change over a variety of system configurations 
and environmental characteristics. Quite simply, the true distribution is neither known nor 
can it be known; yet it can be modeled. To that end, modeling the measurement errors 
with another distribution would prove an interesting test for the SRIF and would be 
relatively straightforward to implement in JBSTAFSim. Investigating the solution 
accuracy and uncertainty using either distributions that are either “heavy tailed” or skewed 
would be worthy endeavors. Likewise, extending JBSTAFSim’s solver to create a non¬ 
diagonal measurement error covariance matrix might prove useful in any of these analyses. 
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Of the many possibilities that exist for the extension of JBSTAFSim, one final 
possibility could be of particular interest. Given its design and network aware capabilities, 
JBSTAFSim’s existing modular design could be modified so that it receives actual TDOA 
and FDOA measurements and sensor information via a network. JBSTAFSim’s solver 
would then compute the geo-location and display the results on a local computer or back 
to the network to be displayed on another host’s computer using only a Web browser. 
While achieving this vision requires the resolution of many security issues. Public Key 
Encryption and secure sockets capabilities within the Java language itself are already being 
implemented in the latest version of Java. Further, this concept agmn shows the “flat” 
architecture that is innate in models created in Java. With very little cost to the user or 
developer, rich capabilities are immediately available to the model which can fimction 
across many platforms in a very robust nature. 
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APPENDIX A. SEQUENTIAL SQUARE ROOT DATA PROCESSING 

ALGORITHM AND THEORY 


The purpose of this appendix is to describe both the underlying theory behind the 
square root information filter (SREF) and the sequential data processing algorithm which 
exploits this theory. A basic understanding of solving linear least squares problems is 
assumed, but is not essential. The primary reference of this section is Reference 16 thus, 
notation will be presented as similarly as possible. In essence this section will provide a 
quick summary reference of the SRIF algorithm as described in detail in Reference 16. 
For the most part, all applicable ideas and examples presented here will be identical to 
those in Reference 16. 

A. CONTEXT OF THE PROBLEM AND REVIEW OF LEAST SQUARES [REF. 
16 , 14 - 16 ] 

With regard to geolocation in a F/TDOA problem, the most basic sequential linear 
least squares system of equations is defined as; 



Where represents the TDOA partials or mixed F/TDOA partials evaluated 

at the state update . This A/xA/ matrix is multiplied by the state step N vector (x- x^,) 

which is approximately equal to the M data vector z minus the TDOA or F/TDOA 
fiinctional A/vector evaluated atX(,. 
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From Chapter FV and for reference, the TDOA/FDOA equations are provided 

below. 


Observation equations; 


c TDOA,^(r) = R^''> -R^'> 

c 


Avf 


Measurement covariance matrix: 


M = 


^ ^TDOA 
0 


0 




T 


<7 


2 

FDOA 


Linearized least squares equation; 


AvO) Avf^ 

j^U) ^(0 J 


0 



-To 


c-TDOA,! y 



-^■FDOA, 


^Av«-Av^V 


Direction cosines are denoted by a, and R is the observer-target separation. In 
the state vector, dx and dv refer to the offset from the position and velocity, respectively, 
of the nominal state, and the subscripts p and 1 refer to components parallel and 
perpendicular, respectively, to the line of sight. The symbols c and fj. refer to the speed 
of light and the frequency of the transmitter. Of course, the dimension of the state here is 
6, not 2 as it appears, because all the elements on the left side of the equation are 
vectors. [Ref 19] 

For simplicity, let us now consider A-Al as the set of equations; 

Ax+ v = z A-A2 


76 

















Here the partials matrix in A-Al is represented by the MxN matrix A ; the state step by the 
vector X; and the functional and data difference by is represented by the vector 2 . We also 
introduce v to represent the M vector of observation errors. The least squares solution to 
minimize the mean square observation error is defined as; 

A-A3 


Or, in terms ofx; 


y(x) = (2 - Ax) ^ (2 - Ax) A-A4 

Thus, J{x) is non-negative and quadratic in the components of x so that the necessary and 
sufiBcient conditions for a minimum require that x satisfy the normal equations: 

A^Ax=A^z A-A5 

A-A5 leads to the well-known least-squares solution x^: 

Xi^={A^Ay'A^z A-A6 

Statistically, if we assume, for now, that the observation errors, v, are random 
variables that have been normalized so that; 

£(v) = 0 and E[v^v) = A-A7 


A-A5 then becomes; 


A^AXi^ = A^z = A^Ax + A^v A-A8 

A-A8 allows us to make two important conclusions about our least squares 
solution. First, that e(x^) = x; thus it x^ is an unbiased estimate. Second, the estimate 

error covariance can be evaluated by: 
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{a'a)e[{x,-x){x^ - xyp'A) = A^e{vv^)A 


A-A9 


A-A9 simplifies so that we get a recognizable estimate for P^, defined as estimate 
error covariance: 




A-AIO 


Where is defined as the information matrix, A^^. 


B. A PRIORI INFORMATION [REF. 16,16-17] 

Now suppose that in addition to the linear system A-A2, there exists an a priori 
unbiased estimate of x and A, defined as x and A. This information can be included in 
the least squares solution by changing A-A4 to: 

= (jc-J^)^A(x-x) + (2-^x)^(2-y4j(:) A-Bl 

Again, like A-A5, the minimizing argument for y,(x), now define as , is defined 
by satisfying the normal equations: 

(a + A^a)x^ = 7^ + A^z A-B2 

From A-B2 and being reminded that (x - x) and v both have mean zero, a similar 
construct to A-A8 can be written to prove that Xj^is unbiased. Further, assuming that the 

observation “noise” and the a priori estimate errors are independent (p[(x-x)v^] = o), 
then like A-A9 and A-AIO we obtain the estimate error covariance: 

P^=£[(.-i,)(r-i,f] = (A + .4-.4)-' A-B3 

Notice how A-B2 and A-B3 reduce to A-A5 and A-AIO when A = 0, i.e. there is 
no a priori information. 
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C. A PRIORI INFORMATION AS ADDITIONAL OBSERVATIONS [REF. 16, 
17-18) 

Given that factorization of a positive semi-definite matrix is always possible, factor 
A so that A = R^R. R^ is called the square root of A. Square root matrices are not 
unique, but are related by orthogonal transformations. For example, if 5', and •S' 2 are 
square roots of a positive semi-definite matrix, then P = S^Sf = , there exists a 

matrix T such that $2 = SyT and TT^ = T^T = I. 

With this factorization in mind, let z = Rx and write; 

(x - x)^ A(x - x) = (z - ^x)^(z - .Rx) A-C1 


Compare these results to J(x) in A-A4 and one can glean that a priori estimate- 

covariance information can be interpreted as additional observations in the least squares 
solution. Therefore, y,(x) in A-Bl can be viewed as simply applying the least squares 

functional J(x) to the augmented system: 



A-C2 


A-C2 is a recursive result since the previous estimate and covariance are combined 
with the new data to form an updated estimate and covariance. Numerically significant is 
that the estimate covariance results do not depend on which square root was used— 
allowing the selection of square root matrices that are simple to compute and enhance 
numerical accuracy. 


D. HOUSEHOLDER ORTHOGONAL TRANSFORMATIONS [REF. 16,57-62] 

In section C, we made a brief reference to orthogonal transformations. SRIF 
depends on the Householder orthogonal transformation. An understanding of this 
transformation is vital in order to understand the SRIF mechanization. 

First, a review of orthogonal matrices. A matrix T is orthogonal if TT^ = I. 
Four properties useful to know about orthogonal transformations are; 
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A. If 7^ and are orthogonal, so is T^T^. 

B. r' = 7’^ and T'^T=I. 

C. For any M vector j, ||7j|| = |[y|| where HI is the Euclidean norm. 

D. If vis an M vector of random variables with 7^(v) = 0 and E{yv^) = I, 
then V = Tv has the same properties as v, namely; ^Cv) = 0 and E{yv^) = / . 

Recall from A-A4 that for MxNA, J{x) ={z- AxY {z - Ax) = \\z - Axf. Define 
rto be a MxM orthogonal matrix and, using Property C from above, we have: 

J{x) = \\T{z-Ax)f = %Tz)-{TA)xf A-Dl 

Since, x^ is independent of T, T can be chosen so that (TM) has a 
computationally attractive form such as: 


n 

m-n 

Where 7? is an upper triangular matrix. If Zz is partitioned. 




A-D2 


A-D3 


then J(x) can be written as: 


J{x) = ||z, - 7?xf + IjzJ' A-D4 

From A-D4, one can see that the minimizing x^^ must satisfy Rx = z^ and that 
j(xjj) = IIZ 2 II . These results have been proved to be less susceptible to errors due to 
computer round-off in References 17 and 18. From above, we see the applicability of 
orthogonal transformations to the least squares problem. Let us now focus on an 
orthogonal transformation that will assist us in the results from section C. 
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Householder transformations are matrix representations corresponding to the 
geometric notion of a reflection. Let « be a nonzero vector and let be the plane 
perpendicular to m . If j is an arbitrary vector. 


y = (y^u)u + v A-D5 

then the reflection of in the plane is equal to: 

>', = -(/«)«+ v A-D6 

where «is a unit vector in the direction of u and v is that part of y that is orthogonal to 
u. See Figure A-Cl. Eliminating v from both equations gives: 

yr=y- 2[y^i^u = (/ - Puu^)y A-D7 


where (3 = r-jj. The matrix T = I - is the elementaiy Householder transformation. 

Il«ll 

Figure A-Cl shows the geometry of the Householder transformation. 



Five properties of this Householder transformation, 7^, are important to us. 

1. \% symmetric. 

2. = I, \s idempotent. (Combined with Property 1 shows that is 
orthogonal.) 
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3. If 1 /( 7 ) = 0 then {Tj^J) = >^ 7 ), rather, if the^h component of u is zero, then 
7^ leaves the 7 th component of y unchanged. 

4. If ul.y then T^y = y. 

5- Tly = y~}^, where y ^ly^uju^u. This is practically a time and storage 
saver since forming and storing 7^ prior to computing T^y requires more 
computation than the direct evaluation using this property. If M is large this 
could be particularly helpful. 

6 . Let <7 = sgn(>{l))-( 7 ^^j)^ and define w(l) =>>(l) + c7, ^( 7 ) = 7 ^ 7),7 > 1. 

Then Tj^ = -oe, and iju^u = l/(aw(l)). 

Property 6 allows us to choose a direction of u so that y^ = lies along e, 
which gives = (cr, 0 ,..., 0 )^. This is the first step in matrix triangularization. Since 
Property C in our orthogonal transformation review assures us that the length of y is 

invariant, therefore It^tI = l<^l = • The direction of u can be obtained from 

Property 5 above: u = const{y - oe^, which explains Property 6 again-= 2 ot/( 1) . 
Refer back to Figure A-Cl; note y,-y is orthogonal to U1 and so is parallel to u . Also 
note, in Property 6 , that the sign of a is chosen to maximize |w(l)| this assists in reducing 
numerical problems when ^ is applied to the columns of A as l/a 7 /(l). 

If Property 6 is applied to an MxN matrix A (with y as its first column) the first 
step of a matrix triangularization results. 


1 n-\ 



A-D8 


Note that via Property 6, s and A are computed directly fi'om A and the matrix 
7^ is implicit—never calculated. Repeated application of Property 6 and the results of 

A-C8 to the columns of A results in a series of orthogonal transformations which 
combine into one orthogonal transformation, T (by Property A of orthogonal 
transformations) and can be written as: 
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A-D9 


TA = 


0 


a, 


a. 


5 . a. 


^i+l 


E. THE SRIF DATA PROCESSING ALGORITHM [REF. 16, 69-72] 

The last four sections now culminate into the SRIF Data Processing Algorithm. 
This method is reputed to be more accurate (less susceptible to the effects of computer 
round-off errors) and stable (accumulated round-off errors will not cause the algorithm to 
diverge) than other conventional non-square root methods. 

Begin by supposing that we have an a priori array from A-Cl [./??] corresponding 

to the equation z = Rx + v where v has zero mean and unity covariance and R is NxN. 
We are interested in constructing the least squares solution to the a priori data equation 
and the new measurements z- Ax+ v. In section C equation A-C2, the least squares 
solution was described as; 



A-El 


In section D, this was seen to be equivalent to finding the least squares solution to: 


T 


R 

A 


x=T 


z 

z 



A-E2 


Further, section D showed that T can be chosen so that; 



A-E3 
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where /? is an upper triangular matrix and we set: 



Z 


f 


V 


V 

T 




and T 





z 


e 


V 




Then A-E2 becomes: 

Rx = z-v 
and 

0 = e-v^ 


A-E4 

A-E5 


A-E4 is in the form of a data equation and can now act as the a priori data 
equation that will be combined with the next set of observations. The e in A-E5 is the 
error in the least squares fit; recall from A-D4: 


J(x) = ||^-z||^ +||e|p 


A-E6 


Thus, l|ef is the sum of squares residual error corresponding to the least squares solution. 
The algorithm then requires the construction of the augmented information 
~R r1 

array and then can be expressed as follows (letting the a priori information array 


equal to 


Tj 


1 

1 

1 

T 


1- 

1- 

> 

1 


-1 

o 

_1 




A-E7 


where [.^o] = [^^o] 7^ orthogonal to triangularize the matrix 

Xj = R'/zj is the least squares estimate corresponding to the a priori information and the 
measurement set, and Pj = R~^R~^ is the covariance of this estimate. Note that because 
Rj is triangular, the calculation of these two estimates can be easily accomplished using 
back substitution methods. 




Thus, 
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F. WHITENING OBSERVATION ERRORS [REF. 16, 47-49] 


Early in section A, equation A-A7 asserted thatE(v^v) = I We now consider 
the case where where P^ is positive definite. We know, then, that we can 

write this matrix as: 


Py = PyJPy A-F1 

where is a lower triangular square root of . Now we can multiply the observations 
by Z,”‘ and the observation set will have unit covariance observation error. Thus the 
familiar Ax+ v = z is transformed to Ax+v =z with: 

2 = Z;'2, A=L~jA, v = LlW 


In the TDOA and TDOA/FDOA case we assume that the observations are uncorrelated, 
but the ith observation error has variance af. Thus, Z,'* = Dia^ / cr,,...1 / <j„). 

Therefore, our augmented matrix takes on the form: 


R z 

L-JA L-Jz 


A-Fl 


As a side note. Reference 16 includes a description of the case where observations are 
correlated. 
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APPENDIX B. MIXED TDOA/FDOA SCALING 


The purpose of this appendix is to discuss the scaling of mixed TDOA and FDOA 
partials, covariance and data in the linear least squares formulation. What follows is a 
summary of and the examples directly from an unpublished work by Dr. Brian Tolman at 
ARL:UT, Reference 19. 

First, consider the usual linear least squares formulation with measurement 
covariance M as described in Ref 19; 


AX = B 

B-1 

X = (A^M-'Ay'A^M-'B 

B-2 

C = (A^M-'Ay' 

B-3 


where. 



B = D-FiX,) 

Dr. Tolman then considers two arbitrary diagonal matrices to re-scale the problem. 


S, = diagonal(mxm) 

B-6 

S 2 = diagonalinxn) 

B-7 

A = S.AS^ 

B-8 

x=syx 

B-9 

B = S,B 

B-10 

M=S,MS, 

B-11 


The problem and its solution then become 
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AX=B 

X= {a^M-'aY'a^M-'B = S-'{A^M-'Ay A^M-'B = S;'X 
C = [A^M-'Ay =S-'{A^M-'AyS-' = 


B-12 

B-13 

B-14 


[Note that if S is diagonal, the effect of multiplying on the left (SA) is to multiply row (i) 
by S(i,i), and the effect of multiplying on the right is to multiply column (j) by S(j,j).] 

The scaling by S^ can be thought of as multiplying each observation equation by a 

separate constant, as when, for example, the units of the TDOA equation are changed 
from seconds to meters. The scaling by can be thought of as changing the units of the 

elements of the state vector (which of course also changes the solution covariance). The 
consistency of the above equations shows that the linear problem may be arbitrarily scaled 
without modifying the solution. In non-linear problems as well, the scaling by clearly 

acts in the same way, effectively multiplying each observation equation by a (perhaps 
different) constant. However, the scaling by cannot work for non-linear problems, 

except in special cases, because there the matrix A is in fact a matrix of partial derivatives 


^ ^X 


B-15 


and the second equality here would hold only if F(X) were linear in X. The TDOA 
equations, however, are a special case (see below), and the state vector can be scaled, 
with one important caveat. [Ref 19] 

Numerical analysis indicates that the ideal situation is to have a matrix A which is 
“equilibrated,” meaning that 

S h Z J ~ B-16 


That is, minimizing the numerical error incurred while solving the least squares equation 
requires that the sum of absolute values along every row and every column has about the 
same magnitude. While this is possible in principle, in practice it is very difficult to find 
scaling matrices which will do this in the general case. [Ref 19] 
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Consider in detail the scaling of the mixed F/TDOA problem as described in Reference 19. 


Observation equations; 


c TDOA^{r) = R^'^ 

Av^'^ 

f-^.FDOA,^{r,v) = -^ 



B-17 

B-18 


Measurement covariance matrix: 


M = 




TDOA 


0 




FDOA 


B-19 


Linearized least squares equation: 


Avp>' 


Avp^ 


rU) pii) 




c5c 

dv. 


c TDOA,^ ^ 

— ■FDOA, 

^fr 


p(') 


O) ^ 


R>'^-R 


^0 


B-20 


Direction cosines are denoted by a, and R is the observer-target separation. In 
the state vector, dx and dv refer to the offset from the position and velocity, respectively, 
of the nominal state, and the subscripts p and J. refer to components parallel and 

perpendicular, respectively, to the line of sight. Note the scaling of these equations. The 
TDOA equation has units (meters) and typically these values are of order 10^ or more. 
The FDOA equation has units (m/s) and typical values are of order 1. Therefore the 
relative scaling is of order 10“’, and may approach much smaller values. This is 
particularly true when the problem begins to approach singularity. Thus there is a 
significant difference of scale within the mixed TDOA/FDOA problem. Also note that the 
TDOA, FDOA and mixed problems, although non-linear, can be scaled in the sense of 
scaling the state vector, as long as all three components of position atid all three 
components of velocity are scaled alike. To see this it is enough to note that changing the 
scale of the three position (velocity) components, uniformly, does not change the physics 
of the problem, it simply changes the units of length (or time or both). [Ref 19] 
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Thus the F/TDOA problem may be scaled in the following way. Choose to scale 
the problem with 4 constants, one each for the three position coordinates together, the 
three velocity components, the TDOA data and the FDOA data. The scaling matrices are 


f r* o'] 




0 F'V 



0 "^ 

vj 


B-21 

B-22 


where each element represents a diagonal block with constant elements, the blocks in 
have dimension 3, and the blocks in S', cover the TDOA data (L) and FDOA data (V). 
Then the least squares equation becomes 






j^U) j^O) 







j-TDOA,j 

f 



— 

FDOA,. 





1 frVIc J 

V 


_ J^U)] / I 


V 


p 

V J 


B-23 


with measurement covariance 


M = 


^ ^ TDOA 

0 


0 

FDOA 


B-24 


Examination of these equations allows one to identify an appropriate value for 
each scaling constant so as to make the terms of order unity, as much as is possible. The 
length scale L should be the scale of the TDOA data, which is of order the typical 
observer-observer separation distance. The other length scale, /, is of order R , the 
observer-target separation. If these length scales are equal the TDOA problem is well 
conditioned and the partials matrix has all terms of order unity. If they are not, the TDOA 
problem is difficult to solve (the geometry is poor), but this is a problem which scaling is 
not going to address anyway. In any case it is probably simplest and best just to choose L 
= / = an average or other compromise of the two length scales. [Ref 19] 
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In the FDOA portion, V apparently needs to be equal to two velocity scales, both 
the scale of the FDOA data (right hand side of the equation), which has the scale of the 
typical observer-observer relative velocity along the target Une-of-sight, and also the scale 
of the typical observer-target relative velocity perpendicular to the line-of-sight (left hand 
side). Again, if these scales are equal the problem is well conditioned, and if they are not, 
there is nothing to do about it. The second velocity scale is arbitrary, since it appears only 
once, in the partials matrix; it is probably best used to cancel the scaling of the direction 
cosines, i.e. v = V. With this scaling in place, clearly every equation is dimensionless, 
and every term has a magnitude of the order of unity (or at least close). Of particular 
importance is the partials matrix, which is nicely “equilibrated”. This should have a direct 
effect on the numerical stability and accuracy of the estimation process. [Ref 19] 
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APPENDIX C. RELATIVE PRECISION 


Relative precision provides a precise estimate of the mean when the sample size of 
the data is unknown. Relative precision requires that we continue to sample our 
simulation results until we are certain, at some level of confidence, a , that the current 
estimate has a given relative error, . Another way of saying this would be that we are 
(l- «)% confident that the percentage error in our estimate, the empirical mean, is 100/ 
percent. From Reference 23, this probabilistic relationship is defined as: 


X-^ 


<r 


= (!-«) 


c-i 


In order to satisfy this relationship, JBSTAFSim samples the miss distance until there have 
been at least 30 observations and the confidence interval half-width of the estimate 
diveded by the empirical mean is less than or equal the percentage error m our estimate. 
Rather, calling the empirical variance and the lOOAth quantile for a standard 

Normal random variable, this terminating condition is defined as: 




n 


X 


r 


\-\-Y 


C-2 


Thus, the left hand side of C-2 represents an estimate of the actual relative error. 
Once y is set, the simulation continues until this estimate satisfies the above inequality 

which means the point estimate X has a relative error of at most //(l->') with 
probability of approximately (l-a). The derivation of C-1 and C-2 can be found in 
Reference 23. 

Finally, it should be noted that our estimate, X, is actually distributed with a 
Student-/ distribution with n degrees of freedom; however since n will be very large 
(greater than 1000), the quantile from the standard Normal distribution is more than 
adequate to calculate the estimate confidence interval half-width. 
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